hair_cortisol_dog_combined

Check working directory

getwd()

## [1] "/Users/alexg/R files/hair_cortisol/skew-normal FINAL"

Load packages

library(readxl)
library(psych)
library(dlookr)

## Registered S3 methods overwritten by 'dlookr':
##   method          from  
##   plot.transform  scales
##   print.transform scales

## 
## Attaching package: 'dlookr'

## The following object is masked from 'package:psych':
## 
##     describe

## The following object is masked from 'package:base':
## 
##     transform

library(vtable)

## Loading required package: kableExtra

library(dplyr)

## 
## Attaching package: 'dplyr'

## The following object is masked from 'package:kableExtra':
## 
##     group_rows

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

library(reshape)

## 
## Attaching package: 'reshape'

## The following object is masked from 'package:dplyr':
## 
##     rename

library(ggplot2)

## 
## Attaching package: 'ggplot2'

## The following objects are masked from 'package:psych':
## 
##     %+%, alpha

library(brms)

## Loading required package: Rcpp

## Loading 'brms' package (version 2.22.0). Useful instructions
## can be found by typing help('brms'). A more detailed introduction
## to the package is available through vignette('brms_overview').

## 
## Attaching package: 'brms'

## The following object is masked from 'package:psych':
## 
##     cs

## The following object is masked from 'package:stats':
## 
##     ar

library(rethinking)

## Loading required package: cmdstanr

## This is cmdstanr version 0.8.0

## - CmdStanR documentation and vignettes: mc-stan.org/cmdstanr

## - CmdStan path: /Users/alexg/.cmdstan/cmdstan-2.36.0

## - CmdStan version: 2.36.0

## Loading required package: posterior

## This is posterior version 1.6.1

## 
## Attaching package: 'posterior'

## The following object is masked from 'package:dlookr':
## 
##     entropy

## The following objects are masked from 'package:stats':
## 
##     mad, sd, var

## The following objects are masked from 'package:base':
## 
##     %in%, match

## Loading required package: parallel

## rethinking (Version 2.42)

## 
## Attaching package: 'rethinking'

## The following objects are masked from 'package:brms':
## 
##     LOO, stancode, WAIC

## The following objects are masked from 'package:psych':
## 
##     logistic, logit, sim

## The following object is masked from 'package:stats':
## 
##     rstudent

library(loo)

## This is loo version 2.8.0

## - Online documentation and vignettes at mc-stan.org/loo

## - As of v2.0.0 loo defaults to 1 core but we recommend using as many as possible. Use the 'cores' argument or set options(mc.cores = NUM_CORES) for an entire session.

## 
## Attaching package: 'loo'

## The following object is masked from 'package:rethinking':
## 
##     compare

library(priorsense)
library(tidyverse)

## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ lubridate 1.9.4     ✔ tibble    3.2.1
## ✔ purrr     1.0.4     ✔ tidyr     1.3.1
## ✔ readr     2.1.5

## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ ggplot2::%+%()      masks psych::%+%()
## ✖ ggplot2::alpha()    masks psych::alpha()
## ✖ tidyr::expand()     masks reshape::expand()
## ✖ tidyr::extract()    masks dlookr::extract()
## ✖ dplyr::filter()     masks stats::filter()
## ✖ dplyr::group_rows() masks kableExtra::group_rows()
## ✖ dplyr::lag()        masks stats::lag()
## ✖ purrr::map()        masks rethinking::map()
## ✖ reshape::rename()   masks dplyr::rename()
## ✖ lubridate::stamp()  masks reshape::stamp()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

library(vioplot)

## Loading required package: sm
## Package 'sm', version 2.2-6.0: type help(sm) for summary information
## 
## Attaching package: 'sm'
## 
## The following object is masked from 'package:dlookr':
## 
##     binning
## 
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## 
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric

library(bayesplot)

## This is bayesplot version 1.12.0
## - Online documentation and vignettes at mc-stan.org/bayesplot
## - bayesplot theme set to bayesplot::theme_default()
##    * Does _not_ affect other ggplot2 plots
##    * See ?bayesplot_theme_set for details on theme setting
## 
## Attaching package: 'bayesplot'
## 
## The following object is masked from 'package:posterior':
## 
##     rhat
## 
## The following object is masked from 'package:brms':
## 
##     rhat

library(bayestestR)

Load data

df <- read_xlsx("hair_cort_dog_all.xlsx", col_types = c("text", "text",  
                               "text", "text", "text", "text",
                               "text", "numeric","text", "skip",
                               "numeric", "skip", "skip", 
                               "numeric", "skip"))
df <- as.data.frame(df)

INITIAL DATA PLOTTING AND EXPLORATION

Check characteristics of df

dim(df) # will tell you how many rows and columns the dataset has

## [1] 73 11

class(df) # will tell you what data structure has the dataset been assigned

## [1] "data.frame"

Explore the dataset to understand its structure.

head(df)

##   number   group visit season breed_group coat_colour    sex age comorbidity
## 1     c1 stopped    v0 winter         ret        dark   Male  43         yes
## 2     c2 stopped    v0 autumn         mix        dark   Male 105         yes
## 3     c3 stopped    v0 spring        ckcs         mix Female 117         yes
## 4     c4 stopped    v0 summer         ret        dark Female 108         yes
## 5     c5 stopped    v0 summer         ret        dark Female 110         yes
## 6     c6 stopped    v0 winter         mix       light Female 120         yes
##   fat_percent cortisol
## 1    52.21393 4.924220
## 2    38.52059 7.304202
## 3    46.94916 1.590000
## 4    44.46813 0.861570
## 5    39.59363 6.217317
## 6          NA 4.426785

1. Get summary stats for numeric data

numeric_df <- Filter(is.numeric, df)
describe(numeric_df) # the describe function in psych provides summary stats

## # A tibble: 3 × 26
##   described_variables     n    na  mean    sd se_mean   IQR skewness kurtosis
##   <chr>               <int> <int> <dbl> <dbl>   <dbl> <dbl>    <dbl>    <dbl>
## 1 age                    73     0 95.8  35.6     4.16 44      -0.104 -0.00589
## 2 fat_percent            55    18 40.5   7.82    1.05  7.82   -0.294  1.12   
## 3 cortisol               73     0  8.11 16.5     1.93  5.43    4.05  18.7    
## # ℹ 17 more variables: p00 <dbl>, p01 <dbl>, p05 <dbl>, p10 <dbl>, p20 <dbl>,
## #   p25 <dbl>, p30 <dbl>, p40 <dbl>, p50 <dbl>, p60 <dbl>, p70 <dbl>,
## #   p75 <dbl>, p80 <dbl>, p90 <dbl>, p95 <dbl>, p99 <dbl>, p100 <dbl>

2. Check normality of all numeric variables

a. graphical assessment

plot_normality(numeric_df)

b. shapiro-wilk test

apply(numeric_df, 2, shapiro.test)

## $age
## 
##  Shapiro-Wilk normality test
## 
## data:  newX[, i]
## W = 0.97361, p-value = 0.1288
## 
## 
## $fat_percent
## 
##  Shapiro-Wilk normality test
## 
## data:  newX[, i]
## W = 0.97956, p-value = 0.4692
## 
## 
## $cortisol
## 
##  Shapiro-Wilk normality test
## 
## data:  newX[, i]
## W = 0.46269, p-value = 6.756e-15

c. repeat Q-Q plots with transformed data

i. log(cortisol)

qqnorm(df$cortisol)
qqline(df$cortisol, col = "red")

qqnorm(log(df$cortisol))
qqline(log(df$cortisol), col = "red")

ii Shapiro test for log cortisol

shapiro.test(log(df$cortisol))

## 
##  Shapiro-Wilk normality test
## 
## data:  log(df$cortisol)
## W = 0.94725, p-value = 0.004126

3. Check data numerically

summary(df$cortisol)

##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
##   0.4141   1.4119   2.3331   8.1089   6.8455 104.6172

summary(log(df$cortisol))

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## -0.8817  0.3449  0.8472  1.1816  1.9236  4.6503

a. log-transform cortisol

df$lgCort <- log(df$cortisol)
summary(df$lgCort)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## -0.8817  0.3449  0.8472  1.1816  1.9236  4.6503

i. visualise

hist(df$lgCort)

b. Create simple category name for breed and convert to factor

df$breed <- df$breed_group
df$breed <- factor(df$breed, levels = c("mix", "ckcs", "pug", "ret", "other"))
head(df$breed)

## [1] ret  mix  ckcs ret  ret  mix 
## Levels: mix ckcs pug ret other

c. Make light hair colour the reference category

df$coat_colour <- factor(df$coat_colour, levels = c("light", "mix", "dark"), ordered = FALSE)
head(df$coat_colour)

## [1] dark  dark  mix   dark  dark  light
## Levels: light mix dark

4. Generate data summary

sumtable(df)

Summary Statistics

Variable	N	Mean	Std. Dev.	Min	Pctl. 25	Pctl. 75	Max
group	73
… completed	42	58%
… stopped	31	42%
visit	73
… v0	52	71%
… v1	21	29%
season	73
… autumn	21	29%
… spring	14	19%
… summer	22	30%
… winter	16	22%
breed_group	73
… ckcs	7	10%
… mix	16	22%
… other	26	36%
… pug	7	10%
… ret	17	23%
coat_colour	73
… light	27	37%
… mix	16	22%
… dark	30	41%
sex	73
… Female	43	59%
… Male	30	41%
age	73	96	36	16	73	117	182
comorbidity	73
… no	15	21%
… yes	58	79%
fat_percent	55	40	7.8	18	37	45	61
cortisol	73	8.1	16	0.41	1.4	6.8	105
lgCort	73	1.2	1.2	-0.88	0.34	1.9	4.7
breed	73
… mix	16	22%
… ckcs	7	10%
… pug	7	10%
… ret	17	23%
… other	26	36%

only stopped data

df_stopped <- df[df$group == "stopped", ]
sumtable(df_stopped, add.median = TRUE)

Summary Statistics

Variable	N	Mean	Std. Dev.	Min	Pctl. 25	Pctl. 50	Pctl. 75	Max
group	31
… stopped	31	100%
visit	31
… v0	31	100%
season	31
… autumn	10	32%
… spring	9	29%
… summer	5	16%
… winter	7	23%
breed_group	31
… ckcs	5	16%
… mix	6	19%
… other	10	32%
… pug	3	10%
… ret	7	23%
coat_colour	31
… light	8	26%
… mix	10	32%
… dark	13	42%
sex	31
… Female	17	55%
… Male	14	45%
age	31	90	34	16	71	98	110	155
comorbidity	31
… no	1	3%
… yes	30	97%
fat_percent	26	41	5.5	31	38	42	44	54
cortisol	31	11	22	0.46	1.8	2.8	7.2	105
lgCort	31	1.4	1.3	-0.77	0.57	1	2	4.7
breed	31
… mix	6	19%
… ckcs	5	16%
… pug	3	10%
… ret	7	23%
… other	10	32%

only completed data

df_completed <- df[df$group == "completed", ]
df_compv0 <- df_completed[df_completed$visit == "v0", ]
df_compv1 <- df_completed[df_completed$visit == "v1", ]
sumtable(df_compv0, add.median = TRUE)

Summary Statistics

Variable	N	Mean	Std. Dev.	Min	Pctl. 25	Pctl. 50	Pctl. 75	Max
group	21
… completed	21	100%
visit	21
… v0	21	100%
season	21
… autumn	3	14%
… spring	4	19%
… summer	8	38%
… winter	6	29%
breed_group	21
… ckcs	1	5%
… mix	5	24%
… other	8	38%
… pug	2	10%
… ret	5	24%
coat_colour	21
… light	10	48%
… mix	3	14%
… dark	8	38%
sex	21
… Female	13	62%
… Male	8	38%
age	21	91	36	29	67	96	109	165
comorbidity	21
… no	7	33%
… yes	14	67%
fat_percent	17	45	6.4	37	40	43	48	61
cortisol	21	7.6	14	0.69	1.7	2.2	6.8	59
lgCort	21	1.1	1.2	-0.38	0.52	0.77	1.9	4.1
breed	21
… mix	5	24%
… ckcs	1	5%
… pug	2	10%
… ret	5	24%
… other	8	38%

sumtable(df_compv1, add.median = TRUE)

Summary Statistics

Variable	N	Mean	Std. Dev.	Min	Pctl. 25	Pctl. 50	Pctl. 75	Max
group	21
… completed	21	100%
visit	21
… v1	21	100%
season	21
… autumn	8	38%
… spring	1	5%
… summer	9	43%
… winter	3	14%
breed_group	21
… ckcs	1	5%
… mix	5	24%
… other	8	38%
… pug	2	10%
… ret	5	24%
coat_colour	21
… light	9	43%
… mix	3	14%
… dark	9	43%
sex	21
… Female	13	62%
… Male	8	38%
age	21	109	35	42	94	105	129	182
comorbidity	21
… no	7	33%
… yes	14	67%
fat_percent	12	33	8.8	18	27	33	38	50
cortisol	21	3.6	3.7	0.41	1.2	1.8	5.6	15
lgCort	21	0.85	0.94	-0.88	0.21	0.58	1.7	2.7
breed	21
… mix	5	24%
… ckcs	1	5%
… pug	2	10%
… ret	5	24%
… other	8	38%

## three way cross tabs (xtabs) and flatten the table
ftable(xtabs(~ visit + group, data = df))

##       group completed stopped
## visit                        
## v0                 21      31
## v1                 21       0

5. Visualise associations

a. Between Cortisol and sex with a violin plot (vioplot package)

par(mfrow = c(1,1))
vioplot(cortisol ~ sex, col = "firebrick",
        data = df)

b. Between lgCortisol and sex with a violin plot (vioplot package)

par(mfrow = c(1,1))
vioplot(lgCort ~ sex, col = "lemonchiffon",
        data = df)

c. Between lgCortisol and breed

i. with a violin plot (vioplot package)

par(mfrow = c(1,1))
vioplot(lgCort ~ breed, col = "firebrick",
        data = df)

ii. with stripchart

stripchart(lgCort ~ breed, vertical = TRUE, method = "jitter",
           col = "steelblue3", data = df, pch = 20)

d. Between Cortisol and comorbidity with a violin plot (vioplot package)

par(mfrow = c(1,1))
vioplot(cortisol ~ comorbidity, col = "steelblue",
        data = df)

e. Between log(cortisol) and comorbidity with a violin plot (vioplot package)

par(mfrow = c(1,1))
vioplot(lgCort ~ comorbidity, col = "steelblue",
        data = df)

f. Between cortisol and age with a dotplot

plot(cortisol ~ age, col = "steelblue3", data = df, pch = 20)

g. Between log(cortisol) and age with a dotplot

plot(lgCort ~ age, col = "steelblue3", data = df, pch = 20)

h. Between cortisol and body fat with a dotplot

plot(cortisol ~ fat_percent, col = "steelblue3", data = df, pch = 20)

h. Between log(cortisol) and body fat with a dotplot

plot(lgCort ~ fat_percent, col = "steelblue3", data = df, pch = 20)

corr.test(df$lgCort, df$fat_percent)

## Call:corr.test(x = df$lgCort, y = df$fat_percent)
## Correlation matrix 
## [1] 0.06
## Sample Size 
## [1] 55
## These are the unadjusted probability values.
##   The probability values  adjusted for multiple tests are in the p.adj object. 
## [1] 0.65
## 
##  To see confidence intervals of the correlations, print with the short=FALSE option

STANDARDISE DATA FOR MODELLING

1. Standardise body fat

df$sFat <- standardize(df$fat_percent)
summary(df$sFat)

##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max.     NA's 
## -2.89165 -0.43568  0.04814  0.00000  0.56366  2.64873       18

a. visualise standardised lgCort

hist(df$sFat)

2. Standardise Age

df$sAge <- standardize(df$age)
summary(df$sAge)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## -2.2444 -0.6422  0.1448  0.0000  0.5945  2.4215

3. Standardise lg(cortisol)

df$slgCort <- standardize(df$lgC)
summary(df$slgCort)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## -1.7079 -0.6925 -0.2768  0.0000  0.6142  2.8713

a. visualise standardised lgCort

hist(df$slgCort, breaks =20, col = "steelblue3", main = "(a) Histogram of observed log hair cortisol", xlab = "Log hair cortisol (standardised)", xlim = c(-2, 3))

b.Recreate using a skew-normal distribution

x <- seq(-2, 4, length.out = 100)
y <- dskew_normal(x, mu = 0, sigma = 1, alpha = 4)
plot(y ~ x, type = "l", main = "(b) Skew-normal distribution with alpha of 4", ylab = "Density", xlab = "Log hair cortisol (standardised)", xlim = c(-2, 3), lwd = 2, col = "steelblue3")

c. Plot the visualisations side by side

par(mfrow = c(1,2))
hist(df$slgCort, breaks =20, col = "steelblue3", main = "(a) Log hair cortisol (observed)", xlab = "Log hair cortisol (standardised)", xlim = c(-2, 3))
x <- seq(-2, 4, length.out = 100)
y <- dskew_normal(x, mu = 0, sigma = 1, alpha = 4)
plot(y ~ x, type = "l", main = "(b) Skew-normal distribution (alpha=4)", ylab = "Density", xlab = "Log hair cortisol (standardised)", xlim = c(-2, 3), lwd = 2, col = "steelblue3")

c. Plot the visualisations with skew-normal overlaid on histogram

par(mfrow = c(1,1))
# prepare histogram data
hist_data <- df$slgCort

# Prepare skew-normal curve dats
x <- seq(-3, 3, length.out = 100)
y <- dskew_normal(x, mu = 0, sigma = 1, alpha = 4)

# Histogram
hist(hist_data, prob = TRUE, col = "steelblue2", breaks = 20,
     ylim = c(0, 0.8), # Adjust y-axis limits to fit the data
     xlim = c(-3, 3),
     main =  "Observed hair cortisol vs skew-normal", 
     xlab = "Log hair cortisol (standardised)",
     ylab = "Density",
     font.lab = 2)

# Add skew-normal curve
lines(x, y, col = "firebrick3", lwd = 3)

# Add normal distribution for comparison
lines(x, dnorm(x, mean = 0, sd = 1), col = "darkblue", lwd = 2, lty = 2)

Create plot to show prior for body fat

x <- seq(-3, 3, length.out = 100)
y <- dnorm(x, mean = 0, sd = 0.5)
plot(y ~ x, type = "l",
     col = "steelblue3", lwd = 3,
      main =  "Prior for the effect of body fat on log hair cortisol", 
     xlab = "Prior probability distribution",
     ylab = "Density",
     font.lab = 2)
abline(v = 0, lwd = 0.5)
abline(v = 0, col = "firebrick3", lty = 2, lwd = 2)

2. create dataset only containing complete data

df2 <- na.omit(df)

MODEL FOR THE EFFECT OF BODY FAT ON HAIR CORTISOL

1. Model code

model <- brm(slgCort ~ sFat + sAge + breed + sex + comorbidity + (1 | visit), family = skew_normal(), data = df2)

2. Check what priors need to be set

default_prior(slgCort ~ sFat + sAge + breed + sex + comorbidity + (1 | visit),
                   family = skew_normal(),
                   data = df2)

##                    prior     class           coef group resp dpar nlpar lb ub
##             normal(0, 4)     alpha                                           
##                   (flat)         b                                           
##                   (flat)         b      breedckcs                            
##                   (flat)         b     breedother                            
##                   (flat)         b       breedpug                            
##                   (flat)         b       breedret                            
##                   (flat)         b comorbidityyes                            
##                   (flat)         b           sAge                            
##                   (flat)         b        sexMale                            
##                   (flat)         b           sFat                            
##  student_t(3, -0.1, 2.5) Intercept                                           
##     student_t(3, 0, 2.5)        sd                                       0   
##     student_t(3, 0, 2.5)        sd                visit                  0   
##     student_t(3, 0, 2.5)        sd      Intercept visit                  0   
##     student_t(3, 0, 2.5)     sigma                                       0   
##        source
##       default
##       default
##  (vectorized)
##  (vectorized)
##  (vectorized)
##  (vectorized)
##  (vectorized)
##  (vectorized)
##  (vectorized)
##  (vectorized)
##       default
##       default
##  (vectorized)
##  (vectorized)
##       default

Published information about associations with hair cortisol

1. body fat effect

1. Log(hair_cort) greater in humans with obesity (logHC = 0.983) cf normal weight (logHC = 0.862) and overweight (0.904) Ref: Obesity (2017) 25, 539-544. doi:10.1002/oby.21733

exp(0.983) # obese

## [1] 2.672462

exp(0.862) # normal

## [1] 2.367892

exp(0.904) # overweight

## [1] 2.469461

Change from normal to obese = 2.37 to 2.67, so only ~10% increase across the range. Prior allowing beta increase of 0.1 per SD increase of BF would work,meaning that we expect change to be between -0.1 and +0.3 Include regularising function to improve efficiency

Change in cortisol when slgCort increases from mean (0) by 0.1 X 1 SD ~10% increase in cortisol for ~10% increase in fat

mlgC <- mean(df2$lgCort)
sdlgC <- sd(df2$lgCort)
exp(mlgC)

## [1] 3.600573

exp(mlgC + 0.1 * sdlgC)

## [1] 4.083976

Change in fat when increases from mean by 1 SD 10% increase in fat

mF <- mean(df2$fat_percent)
sdF <- sd(df2$fat_percent)
mF

## [1] 40.47973

mF + sdF

## [1] 48.30141

Suggests that a 10% increase in fat would be associated with a 10% increase in mean cortisol

Against this evidence in dogs, BCS had negative impact on log hair cortisol… albeit for dogs with BCS 1-6.

NB Bowland found no age effect. They also found a negative effect of BCS on log hair cortisol (beta -0.03). However, BCS ranged from 1-6 (with only 1 BCS 6), so few overweight…. could suggest poor nutrition and health cf obesity. (Bowland JB. Front. Vet. Sci 2020; 7:565346. doi: 10.3389/fvets.2020.565346)

Nonetheless, safer to use a neutral regularising prior.

2. comorbidity effect

1. hair_cort is greater in dogs with atopy vs control Ref: Vet Dermatol 2023 Aug;34(4):339-347. doi: 10.1111/vde.13151.
   hair cortisol before treatment: 0.605 hair cortisol after treatment: 0.25 treatment leads to ~50% reduction

2. hair_cort greater when dogs are stressed (shelter dogs) Ref: Scientific Reports | (2022) 12:5117 | https://doi.org/10.1038/s41598-022-09140-w t0 16.0 +/- 6.8 t1 21.8 +/- 9.4 Average increase = 30-40%

3. hair_cort greater when behavioural stress (competition) and seasonal affect (but variable) Ref: Scientific Reports | 6:19631 | DOI: 10.1038/srep19631 Greater cortisol in January, although only in competing dogs. Cortisol > in competing cf companion and working 40 pg/mg vs 10-20, so approx double.

4. hair_cort greater when behavioural stress (agility) Ref: PLOS ONE | https://doi.org/10.1371/journal.pone.0216000 yes 0.26 (0.17-0.35) not 0.15 (0.10-0.27) Average increase = 50%

Overall, as not exactly the same as effect of various diseases, set a modest (and regularising) prior… 20% increase in cortisol (0.2 SD) for presence of comorbidity

3. sex effect

In humans, males have > hair cortisol cf females (Binz TM. ForensicSciInt 2018; 284:33–8. doi: 10.1016/j.forsciint.2017.12.032) …but effect is opposite in vervet monlkys (Laudenslager ML. Psychoneuroendocrinology 2012; 37:1736–9. doi: 10.1016/j.psyneuen.2012.03.015). … no effect of sex in a previous dog study (Macbeth BJ. Wildl Soc Bull 2012; 36:747–58. doi: 10.1002/wsb.219) … however, this study did find that neutered dogs had decreased hair cortisol. … as all dogs in the study were neutered, this means that an effect of sex is less likely. Bowland et al, female dogs had > hair cortisol than male dogs, but all dogs were intact (Bowland JB. Front. Vet. Sci 2020; 7:565346. doi: 10.3389/fvets.2020.565346)

Further, this effect was lost when accounting for other effects.

Therefore, use a regularising prior but keep it neutral and broad, to allow the effect to be either way.

NB Bowland found no age effect. They also found a negative effect of BCS on log hair cortisol (beta -0.03). However, BCS ranged from 1-6 (with only 1 BCS 6), so few overweight…. could suggest poor nutrition and health cf obesity. (Bowland JB. Front. Vet. Sci 2020; 7:565346. doi: 10.3389/fvets.2020.565346)

breed effect

No published data about effects of breed on hair cortisol, but this is plausible However, unclear as to which breeds will differ and which way. Therefore, use a regularising prior but keep it neutral and broad.

3. Set priors

# Set individual priors
prior_int <- set_prior("normal(0, 0.5)", class = "Intercept")
prior_b <- set_prior("normal(0, 1)", class = "b") # betas for breed
prior_b_sex <- set_prior("normal(0, 1)", class = "b", coef = "sexMale")
prior_b_co <- set_prior("normal(0.25, 1)", class = "b", coef = "comorbidityyes")
prior_b_f <- set_prior("normal(0, 0.5)", class = "b", coef = "sFat")
prior_b_sAge <- set_prior("normal(0, 0.5)", class = "b", coef = "sAge")
prior_sd <- set_prior("normal(0, 1)", class = "sd")
prior_alpha <- set_prior("normal(4, 2)", class = "alpha")

# Combine priors into list
priors <- c(prior_int, prior_b, prior_b_sex, prior_b_co, prior_b_f, prior_b_sAge, prior_sd, prior_alpha)

NB, as we have comorbidity in the nodek, we allowing sigma to vary across comorbidity.

4. Plot priors

a. Prior for intercept

par(mfrow = c(1,1))
x <- seq(-3, 3, length.out = 100)
y <- dnorm(x, mean = 0, sd = 0.5)
plot(y ~ x, type = "l")

b. Prior for sigma

x <- seq(0, 3, length.out = 100)
y <- dexp(x, rate = 1)
plot(y ~ x, type = "l")

b.ii. deciding on alpha for skew normal distribution

Based on distribution of log normal hair cortisol, expect things to be skewed to the right. Try different levels of alpha for skew normal.

x <- seq(-2, 4, length.out = 100)
y <- dskew_normal(x, mu = 0, sigma = 1, alpha = 4)
plot(y ~ x, type = "l")

Given we have a right skew, alpha needs to be positive and 4 appears to reflect the skew seen in the data… keep sd briad (eg 2) broad to allow some flexibility

biii. Prior for visit sd

x <- seq(0, 5, length.out = 100)
y <- dnorm(x, mean = 0, sd = 1)
plot(y ~ x, type = "l")

c. Prior for breed

x <- seq(-3, 3, length.out = 100)
y <- dnorm(x, mean = 0, sd = 1.0)
plot(y ~ x, type = "l")

d. Prior for beta of body fat

x <- seq(-3, 3, length.out = 100)
y <- dnorm(x, mean = 0, sd = 0.5)
plot(y ~ x, type = "l")

e. Prior for beta of comorbidity

x <- seq(-3, 3, length.out = 100)
y <- dnorm(x, mean = 0.25, sd =1)
plot(y ~ x, type = "l", col = "steelblue3", lwd = 3,
     main = "Prior for comorbidity effect on log hair cortisol",
     xlab = "Prior probability distribution",
     ylab = "Density",
     font.lab = 2)
abline(v = 0, lwd = 0.5)
abline(v = 0.25, col = "firebrick3", lty = 2, lwd = 2)

c. Prior for sex

x <- seq(-3, 3, length.out = 100)
y <- dnorm(x, mean = 0, sd = 1)
plot(y ~ x, type = "l")

5. Run model

Increased adapt_delta >0.8 (0.9999 here), as had divergent transitions

set.seed(666)
model <- brm(bf(slgCort ~ sFat + sAge + breed + sex + comorbidity + (1 | visit),
                   sigma ~ comorbidity),
                   family = skew_normal(),
                   prior = priors,
                   data = df2,
                   control=list(adapt_delta=0.999999, stepsize = 0.001, max_treedepth =15),
                   iter = 8000, warmup = 2000,
                   cores = 4,
                   save_pars = save_pars(all =TRUE),
                   sample_prior = TRUE)

## Compiling Stan program...

## Trying to compile a simple C file

## Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
## using C compiler: ‘Apple clang version 17.0.0 (clang-1700.0.13.5)’
## using SDK: ‘MacOSX15.5.sdk’
## clang -arch arm64 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG   -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/Rcpp/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/unsupported"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/src/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppParallel/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/rstan/include" -DEIGEN_NO_DEBUG  -DBOOST_DISABLE_ASSERTS  -DBOOST_PENDING_INTEGER_LOG2_HPP  -DSTAN_THREADS  -DUSE_STANC3 -DSTRICT_R_HEADERS  -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION  -D_HAS_AUTO_PTR_ETC=0  -include '/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp'  -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1   -I/opt/R/arm64/include    -fPIC  -falign-functions=64 -Wall -g -O2  -c foo.c -o foo.o
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22:
## In file included from /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/Eigen/Dense:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/Eigen/Core:19:
## /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:679:10: fatal error: 'cmath' file not found
##   679 | #include <cmath>
##       |          ^~~~~~~
## 1 error generated.
## make: *** [foo.o] Error 1

## Start sampling

## Found more than one class "stanfit" in cache; using the first, from namespace 'rethinking'

## Also defined by 'rstan'

## Found more than one class "stanfit" in cache; using the first, from namespace 'rethinking'

## Also defined by 'rstan'

## Found more than one class "stanfit" in cache; using the first, from namespace 'rethinking'

## Also defined by 'rstan'

## Found more than one class "stanfit" in cache; using the first, from namespace 'rethinking'

## Also defined by 'rstan'

## Warning: There were 3 divergent transitions after warmup. See
## https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
## to find out why this is a problem and how to eliminate them.

## Warning: Examine the pairs() plot to diagnose sampling problems

6. Get summary of model

summary(model)

## Warning: There were 3 divergent transitions after warmup. Increasing
## adapt_delta above 0.999999 may help. See
## http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup

##  Family: skew_normal 
##   Links: mu = identity; sigma = log; alpha = identity 
## Formula: slgCort ~ sFat + sAge + breed + sex + comorbidity + (1 | visit) 
##          sigma ~ comorbidity
##    Data: df2 (Number of observations: 55) 
##   Draws: 4 chains, each with iter = 8000; warmup = 2000; thin = 1;
##          total post-warmup draws = 24000
## 
## Multilevel Hyperparameters:
## ~visit (Number of levels: 2) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.39      0.37     0.01     1.40 1.00    10046    10684
## 
## Regression Coefficients:
##                      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## Intercept               -0.57      0.43    -1.44     0.26 1.00    11140
## sigma_Intercept         -0.94      0.34    -1.52    -0.19 1.00    11614
## sFat                     0.33      0.13     0.04     0.56 1.00    14598
## sAge                     0.06      0.13    -0.21     0.32 1.00    20741
## breedckcs               -0.11      0.44    -0.99     0.73 1.00    19374
## breedpug                -0.17      0.44    -1.03     0.71 1.00    13867
## breedret                -0.23      0.33    -0.88     0.43 1.00    14992
## breedother               0.13      0.34    -0.52     0.83 1.00    13159
## sexMale                  0.10      0.24    -0.38     0.58 1.00    20484
## comorbidityyes           0.80      0.25     0.29     1.28 1.00    16599
## sigma_comorbidityyes     1.09      0.36     0.31     1.73 1.00    12298
##                      Tail_ESS
## Intercept               12627
## sigma_Intercept         12112
## sFat                    14824
## sAge                    17258
## breedckcs               17407
## breedpug                13352
## breedret                14849
## breedother              14221
## sexMale                 16305
## comorbidityyes          14612
## sigma_comorbidityyes    12125
## 
## Further Distributional Parameters:
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## alpha     4.22      1.66     1.19     7.69 1.00    18797    16815
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

7. MCMC diagnostics

a. check distributions and trace plots

color_scheme_set("blue")
plot(model)

Looking for hairy caterpillars

b. try a trank plot as well

mcmc_plot(model, type = 'rank_overlay')

c. trace plot of body_fat only for figure using bayesplot package

i. first create array of posteriors for plotting

posterior <- as.array(model)
dimnames(posterior)

## $iteration
##    [1] "1"    "2"    "3"    "4"    "5"    "6"    "7"    "8"    "9"    "10"  
##   [11] "11"   "12"   "13"   "14"   "15"   "16"   "17"   "18"   "19"   "20"  
##   [21] "21"   "22"   "23"   "24"   "25"   "26"   "27"   "28"   "29"   "30"  
##   [31] "31"   "32"   "33"   "34"   "35"   "36"   "37"   "38"   "39"   "40"  
##   [41] "41"   "42"   "43"   "44"   "45"   "46"   "47"   "48"   "49"   "50"  
##   [51] "51"   "52"   "53"   "54"   "55"   "56"   "57"   "58"   "59"   "60"  
##   [61] "61"   "62"   "63"   "64"   "65"   "66"   "67"   "68"   "69"   "70"  
##   [71] "71"   "72"   "73"   "74"   "75"   "76"   "77"   "78"   "79"   "80"  
##   [81] "81"   "82"   "83"   "84"   "85"   "86"   "87"   "88"   "89"   "90"  
##   [91] "91"   "92"   "93"   "94"   "95"   "96"   "97"   "98"   "99"   "100" 
##  [101] "101"  "102"  "103"  "104"  "105"  "106"  "107"  "108"  "109"  "110" 
##  [111] "111"  "112"  "113"  "114"  "115"  "116"  "117"  "118"  "119"  "120" 
##  [121] "121"  "122"  "123"  "124"  "125"  "126"  "127"  "128"  "129"  "130" 
##  [131] "131"  "132"  "133"  "134"  "135"  "136"  "137"  "138"  "139"  "140" 
##  [141] "141"  "142"  "143"  "144"  "145"  "146"  "147"  "148"  "149"  "150" 
##  [151] "151"  "152"  "153"  "154"  "155"  "156"  "157"  "158"  "159"  "160" 
##  [161] "161"  "162"  "163"  "164"  "165"  "166"  "167"  "168"  "169"  "170" 
##  [171] "171"  "172"  "173"  "174"  "175"  "176"  "177"  "178"  "179"  "180" 
##  [181] "181"  "182"  "183"  "184"  "185"  "186"  "187"  "188"  "189"  "190" 
##  [191] "191"  "192"  "193"  "194"  "195"  "196"  "197"  "198"  "199"  "200" 
##  [201] "201"  "202"  "203"  "204"  "205"  "206"  "207"  "208"  "209"  "210" 
##  [211] "211"  "212"  "213"  "214"  "215"  "216"  "217"  "218"  "219"  "220" 
##  [221] "221"  "222"  "223"  "224"  "225"  "226"  "227"  "228"  "229"  "230" 
##  [231] "231"  "232"  "233"  "234"  "235"  "236"  "237"  "238"  "239"  "240" 
##  [241] "241"  "242"  "243"  "244"  "245"  "246"  "247"  "248"  "249"  "250" 
##  [251] "251"  "252"  "253"  "254"  "255"  "256"  "257"  "258"  "259"  "260" 
##  [261] "261"  "262"  "263"  "264"  "265"  "266"  "267"  "268"  "269"  "270" 
##  [271] "271"  "272"  "273"  "274"  "275"  "276"  "277"  "278"  "279"  "280" 
##  [281] "281"  "282"  "283"  "284"  "285"  "286"  "287"  "288"  "289"  "290" 
##  [291] "291"  "292"  "293"  "294"  "295"  "296"  "297"  "298"  "299"  "300" 
##  [301] "301"  "302"  "303"  "304"  "305"  "306"  "307"  "308"  "309"  "310" 
##  [311] "311"  "312"  "313"  "314"  "315"  "316"  "317"  "318"  "319"  "320" 
##  [321] "321"  "322"  "323"  "324"  "325"  "326"  "327"  "328"  "329"  "330" 
##  [331] "331"  "332"  "333"  "334"  "335"  "336"  "337"  "338"  "339"  "340" 
##  [341] "341"  "342"  "343"  "344"  "345"  "346"  "347"  "348"  "349"  "350" 
##  [351] "351"  "352"  "353"  "354"  "355"  "356"  "357"  "358"  "359"  "360" 
##  [361] "361"  "362"  "363"  "364"  "365"  "366"  "367"  "368"  "369"  "370" 
##  [371] "371"  "372"  "373"  "374"  "375"  "376"  "377"  "378"  "379"  "380" 
##  [381] "381"  "382"  "383"  "384"  "385"  "386"  "387"  "388"  "389"  "390" 
##  [391] "391"  "392"  "393"  "394"  "395"  "396"  "397"  "398"  "399"  "400" 
##  [401] "401"  "402"  "403"  "404"  "405"  "406"  "407"  "408"  "409"  "410" 
##  [411] "411"  "412"  "413"  "414"  "415"  "416"  "417"  "418"  "419"  "420" 
##  [421] "421"  "422"  "423"  "424"  "425"  "426"  "427"  "428"  "429"  "430" 
##  [431] "431"  "432"  "433"  "434"  "435"  "436"  "437"  "438"  "439"  "440" 
##  [441] "441"  "442"  "443"  "444"  "445"  "446"  "447"  "448"  "449"  "450" 
##  [451] "451"  "452"  "453"  "454"  "455"  "456"  "457"  "458"  "459"  "460" 
##  [461] "461"  "462"  "463"  "464"  "465"  "466"  "467"  "468"  "469"  "470" 
##  [471] "471"  "472"  "473"  "474"  "475"  "476"  "477"  "478"  "479"  "480" 
##  [481] "481"  "482"  "483"  "484"  "485"  "486"  "487"  "488"  "489"  "490" 
##  [491] "491"  "492"  "493"  "494"  "495"  "496"  "497"  "498"  "499"  "500" 
##  [501] "501"  "502"  "503"  "504"  "505"  "506"  "507"  "508"  "509"  "510" 
##  [511] "511"  "512"  "513"  "514"  "515"  "516"  "517"  "518"  "519"  "520" 
##  [521] "521"  "522"  "523"  "524"  "525"  "526"  "527"  "528"  "529"  "530" 
##  [531] "531"  "532"  "533"  "534"  "535"  "536"  "537"  "538"  "539"  "540" 
##  [541] "541"  "542"  "543"  "544"  "545"  "546"  "547"  "548"  "549"  "550" 
##  [551] "551"  "552"  "553"  "554"  "555"  "556"  "557"  "558"  "559"  "560" 
##  [561] "561"  "562"  "563"  "564"  "565"  "566"  "567"  "568"  "569"  "570" 
##  [571] "571"  "572"  "573"  "574"  "575"  "576"  "577"  "578"  "579"  "580" 
##  [581] "581"  "582"  "583"  "584"  "585"  "586"  "587"  "588"  "589"  "590" 
##  [591] "591"  "592"  "593"  "594"  "595"  "596"  "597"  "598"  "599"  "600" 
##  [601] "601"  "602"  "603"  "604"  "605"  "606"  "607"  "608"  "609"  "610" 
##  [611] "611"  "612"  "613"  "614"  "615"  "616"  "617"  "618"  "619"  "620" 
##  [621] "621"  "622"  "623"  "624"  "625"  "626"  "627"  "628"  "629"  "630" 
##  [631] "631"  "632"  "633"  "634"  "635"  "636"  "637"  "638"  "639"  "640" 
##  [641] "641"  "642"  "643"  "644"  "645"  "646"  "647"  "648"  "649"  "650" 
##  [651] "651"  "652"  "653"  "654"  "655"  "656"  "657"  "658"  "659"  "660" 
##  [661] "661"  "662"  "663"  "664"  "665"  "666"  "667"  "668"  "669"  "670" 
##  [671] "671"  "672"  "673"  "674"  "675"  "676"  "677"  "678"  "679"  "680" 
##  [681] "681"  "682"  "683"  "684"  "685"  "686"  "687"  "688"  "689"  "690" 
##  [691] "691"  "692"  "693"  "694"  "695"  "696"  "697"  "698"  "699"  "700" 
##  [701] "701"  "702"  "703"  "704"  "705"  "706"  "707"  "708"  "709"  "710" 
##  [711] "711"  "712"  "713"  "714"  "715"  "716"  "717"  "718"  "719"  "720" 
##  [721] "721"  "722"  "723"  "724"  "725"  "726"  "727"  "728"  "729"  "730" 
##  [731] "731"  "732"  "733"  "734"  "735"  "736"  "737"  "738"  "739"  "740" 
##  [741] "741"  "742"  "743"  "744"  "745"  "746"  "747"  "748"  "749"  "750" 
##  [751] "751"  "752"  "753"  "754"  "755"  "756"  "757"  "758"  "759"  "760" 
##  [761] "761"  "762"  "763"  "764"  "765"  "766"  "767"  "768"  "769"  "770" 
##  [771] "771"  "772"  "773"  "774"  "775"  "776"  "777"  "778"  "779"  "780" 
##  [781] "781"  "782"  "783"  "784"  "785"  "786"  "787"  "788"  "789"  "790" 
##  [791] "791"  "792"  "793"  "794"  "795"  "796"  "797"  "798"  "799"  "800" 
##  [801] "801"  "802"  "803"  "804"  "805"  "806"  "807"  "808"  "809"  "810" 
##  [811] "811"  "812"  "813"  "814"  "815"  "816"  "817"  "818"  "819"  "820" 
##  [821] "821"  "822"  "823"  "824"  "825"  "826"  "827"  "828"  "829"  "830" 
##  [831] "831"  "832"  "833"  "834"  "835"  "836"  "837"  "838"  "839"  "840" 
##  [841] "841"  "842"  "843"  "844"  "845"  "846"  "847"  "848"  "849"  "850" 
##  [851] "851"  "852"  "853"  "854"  "855"  "856"  "857"  "858"  "859"  "860" 
##  [861] "861"  "862"  "863"  "864"  "865"  "866"  "867"  "868"  "869"  "870" 
##  [871] "871"  "872"  "873"  "874"  "875"  "876"  "877"  "878"  "879"  "880" 
##  [881] "881"  "882"  "883"  "884"  "885"  "886"  "887"  "888"  "889"  "890" 
##  [891] "891"  "892"  "893"  "894"  "895"  "896"  "897"  "898"  "899"  "900" 
##  [901] "901"  "902"  "903"  "904"  "905"  "906"  "907"  "908"  "909"  "910" 
##  [911] "911"  "912"  "913"  "914"  "915"  "916"  "917"  "918"  "919"  "920" 
##  [921] "921"  "922"  "923"  "924"  "925"  "926"  "927"  "928"  "929"  "930" 
##  [931] "931"  "932"  "933"  "934"  "935"  "936"  "937"  "938"  "939"  "940" 
##  [941] "941"  "942"  "943"  "944"  "945"  "946"  "947"  "948"  "949"  "950" 
##  [951] "951"  "952"  "953"  "954"  "955"  "956"  "957"  "958"  "959"  "960" 
##  [961] "961"  "962"  "963"  "964"  "965"  "966"  "967"  "968"  "969"  "970" 
##  [971] "971"  "972"  "973"  "974"  "975"  "976"  "977"  "978"  "979"  "980" 
##  [981] "981"  "982"  "983"  "984"  "985"  "986"  "987"  "988"  "989"  "990" 
##  [991] "991"  "992"  "993"  "994"  "995"  "996"  "997"  "998"  "999"  "1000"
## [1001] "1001" "1002" "1003" "1004" "1005" "1006" "1007" "1008" "1009" "1010"
## [1011] "1011" "1012" "1013" "1014" "1015" "1016" "1017" "1018" "1019" "1020"
## [1021] "1021" "1022" "1023" "1024" "1025" "1026" "1027" "1028" "1029" "1030"
## [1031] "1031" "1032" "1033" "1034" "1035" "1036" "1037" "1038" "1039" "1040"
## [1041] "1041" "1042" "1043" "1044" "1045" "1046" "1047" "1048" "1049" "1050"
## [1051] "1051" "1052" "1053" "1054" "1055" "1056" "1057" "1058" "1059" "1060"
## [1061] "1061" "1062" "1063" "1064" "1065" "1066" "1067" "1068" "1069" "1070"
## [1071] "1071" "1072" "1073" "1074" "1075" "1076" "1077" "1078" "1079" "1080"
## [1081] "1081" "1082" "1083" "1084" "1085" "1086" "1087" "1088" "1089" "1090"
## [1091] "1091" "1092" "1093" "1094" "1095" "1096" "1097" "1098" "1099" "1100"
## [1101] "1101" "1102" "1103" "1104" "1105" "1106" "1107" "1108" "1109" "1110"
## [1111] "1111" "1112" "1113" "1114" "1115" "1116" "1117" "1118" "1119" "1120"
## [1121] "1121" "1122" "1123" "1124" "1125" "1126" "1127" "1128" "1129" "1130"
## [1131] "1131" "1132" "1133" "1134" "1135" "1136" "1137" "1138" "1139" "1140"
## [1141] "1141" "1142" "1143" "1144" "1145" "1146" "1147" "1148" "1149" "1150"
## [1151] "1151" "1152" "1153" "1154" "1155" "1156" "1157" "1158" "1159" "1160"
## [1161] "1161" "1162" "1163" "1164" "1165" "1166" "1167" "1168" "1169" "1170"
## [1171] "1171" "1172" "1173" "1174" "1175" "1176" "1177" "1178" "1179" "1180"
## [1181] "1181" "1182" "1183" "1184" "1185" "1186" "1187" "1188" "1189" "1190"
## [1191] "1191" "1192" "1193" "1194" "1195" "1196" "1197" "1198" "1199" "1200"
## [1201] "1201" "1202" "1203" "1204" "1205" "1206" "1207" "1208" "1209" "1210"
## [1211] "1211" "1212" "1213" "1214" "1215" "1216" "1217" "1218" "1219" "1220"
## [1221] "1221" "1222" "1223" "1224" "1225" "1226" "1227" "1228" "1229" "1230"
## [1231] "1231" "1232" "1233" "1234" "1235" "1236" "1237" "1238" "1239" "1240"
## [1241] "1241" "1242" "1243" "1244" "1245" "1246" "1247" "1248" "1249" "1250"
## [1251] "1251" "1252" "1253" "1254" "1255" "1256" "1257" "1258" "1259" "1260"
## [1261] "1261" "1262" "1263" "1264" "1265" "1266" "1267" "1268" "1269" "1270"
## [1271] "1271" "1272" "1273" "1274" "1275" "1276" "1277" "1278" "1279" "1280"
## [1281] "1281" "1282" "1283" "1284" "1285" "1286" "1287" "1288" "1289" "1290"
## [1291] "1291" "1292" "1293" "1294" "1295" "1296" "1297" "1298" "1299" "1300"
## [1301] "1301" "1302" "1303" "1304" "1305" "1306" "1307" "1308" "1309" "1310"
## [1311] "1311" "1312" "1313" "1314" "1315" "1316" "1317" "1318" "1319" "1320"
## [1321] "1321" "1322" "1323" "1324" "1325" "1326" "1327" "1328" "1329" "1330"
## [1331] "1331" "1332" "1333" "1334" "1335" "1336" "1337" "1338" "1339" "1340"
## [1341] "1341" "1342" "1343" "1344" "1345" "1346" "1347" "1348" "1349" "1350"
## [1351] "1351" "1352" "1353" "1354" "1355" "1356" "1357" "1358" "1359" "1360"
## [1361] "1361" "1362" "1363" "1364" "1365" "1366" "1367" "1368" "1369" "1370"
## [1371] "1371" "1372" "1373" "1374" "1375" "1376" "1377" "1378" "1379" "1380"
## [1381] "1381" "1382" "1383" "1384" "1385" "1386" "1387" "1388" "1389" "1390"
## [1391] "1391" "1392" "1393" "1394" "1395" "1396" "1397" "1398" "1399" "1400"
## [1401] "1401" "1402" "1403" "1404" "1405" "1406" "1407" "1408" "1409" "1410"
## [1411] "1411" "1412" "1413" "1414" "1415" "1416" "1417" "1418" "1419" "1420"
## [1421] "1421" "1422" "1423" "1424" "1425" "1426" "1427" "1428" "1429" "1430"
## [1431] "1431" "1432" "1433" "1434" "1435" "1436" "1437" "1438" "1439" "1440"
## [1441] "1441" "1442" "1443" "1444" "1445" "1446" "1447" "1448" "1449" "1450"
## [1451] "1451" "1452" "1453" "1454" "1455" "1456" "1457" "1458" "1459" "1460"
## [1461] "1461" "1462" "1463" "1464" "1465" "1466" "1467" "1468" "1469" "1470"
## [1471] "1471" "1472" "1473" "1474" "1475" "1476" "1477" "1478" "1479" "1480"
## [1481] "1481" "1482" "1483" "1484" "1485" "1486" "1487" "1488" "1489" "1490"
## [1491] "1491" "1492" "1493" "1494" "1495" "1496" "1497" "1498" "1499" "1500"
## [1501] "1501" "1502" "1503" "1504" "1505" "1506" "1507" "1508" "1509" "1510"
## [1511] "1511" "1512" "1513" "1514" "1515" "1516" "1517" "1518" "1519" "1520"
## [1521] "1521" "1522" "1523" "1524" "1525" "1526" "1527" "1528" "1529" "1530"
## [1531] "1531" "1532" "1533" "1534" "1535" "1536" "1537" "1538" "1539" "1540"
## [1541] "1541" "1542" "1543" "1544" "1545" "1546" "1547" "1548" "1549" "1550"
## [1551] "1551" "1552" "1553" "1554" "1555" "1556" "1557" "1558" "1559" "1560"
## [1561] "1561" "1562" "1563" "1564" "1565" "1566" "1567" "1568" "1569" "1570"
## [1571] "1571" "1572" "1573" "1574" "1575" "1576" "1577" "1578" "1579" "1580"
## [1581] "1581" "1582" "1583" "1584" "1585" "1586" "1587" "1588" "1589" "1590"
## [1591] "1591" "1592" "1593" "1594" "1595" "1596" "1597" "1598" "1599" "1600"
## [1601] "1601" "1602" "1603" "1604" "1605" "1606" "1607" "1608" "1609" "1610"
## [1611] "1611" "1612" "1613" "1614" "1615" "1616" "1617" "1618" "1619" "1620"
## [1621] "1621" "1622" "1623" "1624" "1625" "1626" "1627" "1628" "1629" "1630"
## [1631] "1631" "1632" "1633" "1634" "1635" "1636" "1637" "1638" "1639" "1640"
## [1641] "1641" "1642" "1643" "1644" "1645" "1646" "1647" "1648" "1649" "1650"
## [1651] "1651" "1652" "1653" "1654" "1655" "1656" "1657" "1658" "1659" "1660"
## [1661] "1661" "1662" "1663" "1664" "1665" "1666" "1667" "1668" "1669" "1670"
## [1671] "1671" "1672" "1673" "1674" "1675" "1676" "1677" "1678" "1679" "1680"
## [1681] "1681" "1682" "1683" "1684" "1685" "1686" "1687" "1688" "1689" "1690"
## [1691] "1691" "1692" "1693" "1694" "1695" "1696" "1697" "1698" "1699" "1700"
## [1701] "1701" "1702" "1703" "1704" "1705" "1706" "1707" "1708" "1709" "1710"
## [1711] "1711" "1712" "1713" "1714" "1715" "1716" "1717" "1718" "1719" "1720"
## [1721] "1721" "1722" "1723" "1724" "1725" "1726" "1727" "1728" "1729" "1730"
## [1731] "1731" "1732" "1733" "1734" "1735" "1736" "1737" "1738" "1739" "1740"
## [1741] "1741" "1742" "1743" "1744" "1745" "1746" "1747" "1748" "1749" "1750"
## [1751] "1751" "1752" "1753" "1754" "1755" "1756" "1757" "1758" "1759" "1760"
## [1761] "1761" "1762" "1763" "1764" "1765" "1766" "1767" "1768" "1769" "1770"
## [1771] "1771" "1772" "1773" "1774" "1775" "1776" "1777" "1778" "1779" "1780"
## [1781] "1781" "1782" "1783" "1784" "1785" "1786" "1787" "1788" "1789" "1790"
## [1791] "1791" "1792" "1793" "1794" "1795" "1796" "1797" "1798" "1799" "1800"
## [1801] "1801" "1802" "1803" "1804" "1805" "1806" "1807" "1808" "1809" "1810"
## [1811] "1811" "1812" "1813" "1814" "1815" "1816" "1817" "1818" "1819" "1820"
## [1821] "1821" "1822" "1823" "1824" "1825" "1826" "1827" "1828" "1829" "1830"
## [1831] "1831" "1832" "1833" "1834" "1835" "1836" "1837" "1838" "1839" "1840"
## [1841] "1841" "1842" "1843" "1844" "1845" "1846" "1847" "1848" "1849" "1850"
## [1851] "1851" "1852" "1853" "1854" "1855" "1856" "1857" "1858" "1859" "1860"
## [1861] "1861" "1862" "1863" "1864" "1865" "1866" "1867" "1868" "1869" "1870"
## [1871] "1871" "1872" "1873" "1874" "1875" "1876" "1877" "1878" "1879" "1880"
## [1881] "1881" "1882" "1883" "1884" "1885" "1886" "1887" "1888" "1889" "1890"
## [1891] "1891" "1892" "1893" "1894" "1895" "1896" "1897" "1898" "1899" "1900"
## [1901] "1901" "1902" "1903" "1904" "1905" "1906" "1907" "1908" "1909" "1910"
## [1911] "1911" "1912" "1913" "1914" "1915" "1916" "1917" "1918" "1919" "1920"
## [1921] "1921" "1922" "1923" "1924" "1925" "1926" "1927" "1928" "1929" "1930"
## [1931] "1931" "1932" "1933" "1934" "1935" "1936" "1937" "1938" "1939" "1940"
## [1941] "1941" "1942" "1943" "1944" "1945" "1946" "1947" "1948" "1949" "1950"
## [1951] "1951" "1952" "1953" "1954" "1955" "1956" "1957" "1958" "1959" "1960"
## [1961] "1961" "1962" "1963" "1964" "1965" "1966" "1967" "1968" "1969" "1970"
## [1971] "1971" "1972" "1973" "1974" "1975" "1976" "1977" "1978" "1979" "1980"
## [1981] "1981" "1982" "1983" "1984" "1985" "1986" "1987" "1988" "1989" "1990"
## [1991] "1991" "1992" "1993" "1994" "1995" "1996" "1997" "1998" "1999" "2000"
## [2001] "2001" "2002" "2003" "2004" "2005" "2006" "2007" "2008" "2009" "2010"
## [2011] "2011" "2012" "2013" "2014" "2015" "2016" "2017" "2018" "2019" "2020"
## [2021] "2021" "2022" "2023" "2024" "2025" "2026" "2027" "2028" "2029" "2030"
## [2031] "2031" "2032" "2033" "2034" "2035" "2036" "2037" "2038" "2039" "2040"
## [2041] "2041" "2042" "2043" "2044" "2045" "2046" "2047" "2048" "2049" "2050"
## [2051] "2051" "2052" "2053" "2054" "2055" "2056" "2057" "2058" "2059" "2060"
## [2061] "2061" "2062" "2063" "2064" "2065" "2066" "2067" "2068" "2069" "2070"
## [2071] "2071" "2072" "2073" "2074" "2075" "2076" "2077" "2078" "2079" "2080"
## [2081] "2081" "2082" "2083" "2084" "2085" "2086" "2087" "2088" "2089" "2090"
## [2091] "2091" "2092" "2093" "2094" "2095" "2096" "2097" "2098" "2099" "2100"
## [2101] "2101" "2102" "2103" "2104" "2105" "2106" "2107" "2108" "2109" "2110"
## [2111] "2111" "2112" "2113" "2114" "2115" "2116" "2117" "2118" "2119" "2120"
## [2121] "2121" "2122" "2123" "2124" "2125" "2126" "2127" "2128" "2129" "2130"
## [2131] "2131" "2132" "2133" "2134" "2135" "2136" "2137" "2138" "2139" "2140"
## [2141] "2141" "2142" "2143" "2144" "2145" "2146" "2147" "2148" "2149" "2150"
## [2151] "2151" "2152" "2153" "2154" "2155" "2156" "2157" "2158" "2159" "2160"
## [2161] "2161" "2162" "2163" "2164" "2165" "2166" "2167" "2168" "2169" "2170"
## [2171] "2171" "2172" "2173" "2174" "2175" "2176" "2177" "2178" "2179" "2180"
## [2181] "2181" "2182" "2183" "2184" "2185" "2186" "2187" "2188" "2189" "2190"
## [2191] "2191" "2192" "2193" "2194" "2195" "2196" "2197" "2198" "2199" "2200"
## [2201] "2201" "2202" "2203" "2204" "2205" "2206" "2207" "2208" "2209" "2210"
## [2211] "2211" "2212" "2213" "2214" "2215" "2216" "2217" "2218" "2219" "2220"
## [2221] "2221" "2222" "2223" "2224" "2225" "2226" "2227" "2228" "2229" "2230"
## [2231] "2231" "2232" "2233" "2234" "2235" "2236" "2237" "2238" "2239" "2240"
## [2241] "2241" "2242" "2243" "2244" "2245" "2246" "2247" "2248" "2249" "2250"
## [2251] "2251" "2252" "2253" "2254" "2255" "2256" "2257" "2258" "2259" "2260"
## [2261] "2261" "2262" "2263" "2264" "2265" "2266" "2267" "2268" "2269" "2270"
## [2271] "2271" "2272" "2273" "2274" "2275" "2276" "2277" "2278" "2279" "2280"
## [2281] "2281" "2282" "2283" "2284" "2285" "2286" "2287" "2288" "2289" "2290"
## [2291] "2291" "2292" "2293" "2294" "2295" "2296" "2297" "2298" "2299" "2300"
## [2301] "2301" "2302" "2303" "2304" "2305" "2306" "2307" "2308" "2309" "2310"
## [2311] "2311" "2312" "2313" "2314" "2315" "2316" "2317" "2318" "2319" "2320"
## [2321] "2321" "2322" "2323" "2324" "2325" "2326" "2327" "2328" "2329" "2330"
## [2331] "2331" "2332" "2333" "2334" "2335" "2336" "2337" "2338" "2339" "2340"
## [2341] "2341" "2342" "2343" "2344" "2345" "2346" "2347" "2348" "2349" "2350"
## [2351] "2351" "2352" "2353" "2354" "2355" "2356" "2357" "2358" "2359" "2360"
## [2361] "2361" "2362" "2363" "2364" "2365" "2366" "2367" "2368" "2369" "2370"
## [2371] "2371" "2372" "2373" "2374" "2375" "2376" "2377" "2378" "2379" "2380"
## [2381] "2381" "2382" "2383" "2384" "2385" "2386" "2387" "2388" "2389" "2390"
## [2391] "2391" "2392" "2393" "2394" "2395" "2396" "2397" "2398" "2399" "2400"
## [2401] "2401" "2402" "2403" "2404" "2405" "2406" "2407" "2408" "2409" "2410"
## [2411] "2411" "2412" "2413" "2414" "2415" "2416" "2417" "2418" "2419" "2420"
## [2421] "2421" "2422" "2423" "2424" "2425" "2426" "2427" "2428" "2429" "2430"
## [2431] "2431" "2432" "2433" "2434" "2435" "2436" "2437" "2438" "2439" "2440"
## [2441] "2441" "2442" "2443" "2444" "2445" "2446" "2447" "2448" "2449" "2450"
## [2451] "2451" "2452" "2453" "2454" "2455" "2456" "2457" "2458" "2459" "2460"
## [2461] "2461" "2462" "2463" "2464" "2465" "2466" "2467" "2468" "2469" "2470"
## [2471] "2471" "2472" "2473" "2474" "2475" "2476" "2477" "2478" "2479" "2480"
## [2481] "2481" "2482" "2483" "2484" "2485" "2486" "2487" "2488" "2489" "2490"
## [2491] "2491" "2492" "2493" "2494" "2495" "2496" "2497" "2498" "2499" "2500"
## [2501] "2501" "2502" "2503" "2504" "2505" "2506" "2507" "2508" "2509" "2510"
## [2511] "2511" "2512" "2513" "2514" "2515" "2516" "2517" "2518" "2519" "2520"
## [2521] "2521" "2522" "2523" "2524" "2525" "2526" "2527" "2528" "2529" "2530"
## [2531] "2531" "2532" "2533" "2534" "2535" "2536" "2537" "2538" "2539" "2540"
## [2541] "2541" "2542" "2543" "2544" "2545" "2546" "2547" "2548" "2549" "2550"
## [2551] "2551" "2552" "2553" "2554" "2555" "2556" "2557" "2558" "2559" "2560"
## [2561] "2561" "2562" "2563" "2564" "2565" "2566" "2567" "2568" "2569" "2570"
## [2571] "2571" "2572" "2573" "2574" "2575" "2576" "2577" "2578" "2579" "2580"
## [2581] "2581" "2582" "2583" "2584" "2585" "2586" "2587" "2588" "2589" "2590"
## [2591] "2591" "2592" "2593" "2594" "2595" "2596" "2597" "2598" "2599" "2600"
## [2601] "2601" "2602" "2603" "2604" "2605" "2606" "2607" "2608" "2609" "2610"
## [2611] "2611" "2612" "2613" "2614" "2615" "2616" "2617" "2618" "2619" "2620"
## [2621] "2621" "2622" "2623" "2624" "2625" "2626" "2627" "2628" "2629" "2630"
## [2631] "2631" "2632" "2633" "2634" "2635" "2636" "2637" "2638" "2639" "2640"
## [2641] "2641" "2642" "2643" "2644" "2645" "2646" "2647" "2648" "2649" "2650"
## [2651] "2651" "2652" "2653" "2654" "2655" "2656" "2657" "2658" "2659" "2660"
## [2661] "2661" "2662" "2663" "2664" "2665" "2666" "2667" "2668" "2669" "2670"
## [2671] "2671" "2672" "2673" "2674" "2675" "2676" "2677" "2678" "2679" "2680"
## [2681] "2681" "2682" "2683" "2684" "2685" "2686" "2687" "2688" "2689" "2690"
## [2691] "2691" "2692" "2693" "2694" "2695" "2696" "2697" "2698" "2699" "2700"
## [2701] "2701" "2702" "2703" "2704" "2705" "2706" "2707" "2708" "2709" "2710"
## [2711] "2711" "2712" "2713" "2714" "2715" "2716" "2717" "2718" "2719" "2720"
## [2721] "2721" "2722" "2723" "2724" "2725" "2726" "2727" "2728" "2729" "2730"
## [2731] "2731" "2732" "2733" "2734" "2735" "2736" "2737" "2738" "2739" "2740"
## [2741] "2741" "2742" "2743" "2744" "2745" "2746" "2747" "2748" "2749" "2750"
## [2751] "2751" "2752" "2753" "2754" "2755" "2756" "2757" "2758" "2759" "2760"
## [2761] "2761" "2762" "2763" "2764" "2765" "2766" "2767" "2768" "2769" "2770"
## [2771] "2771" "2772" "2773" "2774" "2775" "2776" "2777" "2778" "2779" "2780"
## [2781] "2781" "2782" "2783" "2784" "2785" "2786" "2787" "2788" "2789" "2790"
## [2791] "2791" "2792" "2793" "2794" "2795" "2796" "2797" "2798" "2799" "2800"
## [2801] "2801" "2802" "2803" "2804" "2805" "2806" "2807" "2808" "2809" "2810"
## [2811] "2811" "2812" "2813" "2814" "2815" "2816" "2817" "2818" "2819" "2820"
## [2821] "2821" "2822" "2823" "2824" "2825" "2826" "2827" "2828" "2829" "2830"
## [2831] "2831" "2832" "2833" "2834" "2835" "2836" "2837" "2838" "2839" "2840"
## [2841] "2841" "2842" "2843" "2844" "2845" "2846" "2847" "2848" "2849" "2850"
## [2851] "2851" "2852" "2853" "2854" "2855" "2856" "2857" "2858" "2859" "2860"
## [2861] "2861" "2862" "2863" "2864" "2865" "2866" "2867" "2868" "2869" "2870"
## [2871] "2871" "2872" "2873" "2874" "2875" "2876" "2877" "2878" "2879" "2880"
## [2881] "2881" "2882" "2883" "2884" "2885" "2886" "2887" "2888" "2889" "2890"
## [2891] "2891" "2892" "2893" "2894" "2895" "2896" "2897" "2898" "2899" "2900"
## [2901] "2901" "2902" "2903" "2904" "2905" "2906" "2907" "2908" "2909" "2910"
## [2911] "2911" "2912" "2913" "2914" "2915" "2916" "2917" "2918" "2919" "2920"
## [2921] "2921" "2922" "2923" "2924" "2925" "2926" "2927" "2928" "2929" "2930"
## [2931] "2931" "2932" "2933" "2934" "2935" "2936" "2937" "2938" "2939" "2940"
## [2941] "2941" "2942" "2943" "2944" "2945" "2946" "2947" "2948" "2949" "2950"
## [2951] "2951" "2952" "2953" "2954" "2955" "2956" "2957" "2958" "2959" "2960"
## [2961] "2961" "2962" "2963" "2964" "2965" "2966" "2967" "2968" "2969" "2970"
## [2971] "2971" "2972" "2973" "2974" "2975" "2976" "2977" "2978" "2979" "2980"
## [2981] "2981" "2982" "2983" "2984" "2985" "2986" "2987" "2988" "2989" "2990"
## [2991] "2991" "2992" "2993" "2994" "2995" "2996" "2997" "2998" "2999" "3000"
## [3001] "3001" "3002" "3003" "3004" "3005" "3006" "3007" "3008" "3009" "3010"
## [3011] "3011" "3012" "3013" "3014" "3015" "3016" "3017" "3018" "3019" "3020"
## [3021] "3021" "3022" "3023" "3024" "3025" "3026" "3027" "3028" "3029" "3030"
## [3031] "3031" "3032" "3033" "3034" "3035" "3036" "3037" "3038" "3039" "3040"
## [3041] "3041" "3042" "3043" "3044" "3045" "3046" "3047" "3048" "3049" "3050"
## [3051] "3051" "3052" "3053" "3054" "3055" "3056" "3057" "3058" "3059" "3060"
## [3061] "3061" "3062" "3063" "3064" "3065" "3066" "3067" "3068" "3069" "3070"
## [3071] "3071" "3072" "3073" "3074" "3075" "3076" "3077" "3078" "3079" "3080"
## [3081] "3081" "3082" "3083" "3084" "3085" "3086" "3087" "3088" "3089" "3090"
## [3091] "3091" "3092" "3093" "3094" "3095" "3096" "3097" "3098" "3099" "3100"
## [3101] "3101" "3102" "3103" "3104" "3105" "3106" "3107" "3108" "3109" "3110"
## [3111] "3111" "3112" "3113" "3114" "3115" "3116" "3117" "3118" "3119" "3120"
## [3121] "3121" "3122" "3123" "3124" "3125" "3126" "3127" "3128" "3129" "3130"
## [3131] "3131" "3132" "3133" "3134" "3135" "3136" "3137" "3138" "3139" "3140"
## [3141] "3141" "3142" "3143" "3144" "3145" "3146" "3147" "3148" "3149" "3150"
## [3151] "3151" "3152" "3153" "3154" "3155" "3156" "3157" "3158" "3159" "3160"
## [3161] "3161" "3162" "3163" "3164" "3165" "3166" "3167" "3168" "3169" "3170"
## [3171] "3171" "3172" "3173" "3174" "3175" "3176" "3177" "3178" "3179" "3180"
## [3181] "3181" "3182" "3183" "3184" "3185" "3186" "3187" "3188" "3189" "3190"
## [3191] "3191" "3192" "3193" "3194" "3195" "3196" "3197" "3198" "3199" "3200"
## [3201] "3201" "3202" "3203" "3204" "3205" "3206" "3207" "3208" "3209" "3210"
## [3211] "3211" "3212" "3213" "3214" "3215" "3216" "3217" "3218" "3219" "3220"
## [3221] "3221" "3222" "3223" "3224" "3225" "3226" "3227" "3228" "3229" "3230"
## [3231] "3231" "3232" "3233" "3234" "3235" "3236" "3237" "3238" "3239" "3240"
## [3241] "3241" "3242" "3243" "3244" "3245" "3246" "3247" "3248" "3249" "3250"
## [3251] "3251" "3252" "3253" "3254" "3255" "3256" "3257" "3258" "3259" "3260"
## [3261] "3261" "3262" "3263" "3264" "3265" "3266" "3267" "3268" "3269" "3270"
## [3271] "3271" "3272" "3273" "3274" "3275" "3276" "3277" "3278" "3279" "3280"
## [3281] "3281" "3282" "3283" "3284" "3285" "3286" "3287" "3288" "3289" "3290"
## [3291] "3291" "3292" "3293" "3294" "3295" "3296" "3297" "3298" "3299" "3300"
## [3301] "3301" "3302" "3303" "3304" "3305" "3306" "3307" "3308" "3309" "3310"
## [3311] "3311" "3312" "3313" "3314" "3315" "3316" "3317" "3318" "3319" "3320"
## [3321] "3321" "3322" "3323" "3324" "3325" "3326" "3327" "3328" "3329" "3330"
## [3331] "3331" "3332" "3333" "3334" "3335" "3336" "3337" "3338" "3339" "3340"
## [3341] "3341" "3342" "3343" "3344" "3345" "3346" "3347" "3348" "3349" "3350"
## [3351] "3351" "3352" "3353" "3354" "3355" "3356" "3357" "3358" "3359" "3360"
## [3361] "3361" "3362" "3363" "3364" "3365" "3366" "3367" "3368" "3369" "3370"
## [3371] "3371" "3372" "3373" "3374" "3375" "3376" "3377" "3378" "3379" "3380"
## [3381] "3381" "3382" "3383" "3384" "3385" "3386" "3387" "3388" "3389" "3390"
## [3391] "3391" "3392" "3393" "3394" "3395" "3396" "3397" "3398" "3399" "3400"
## [3401] "3401" "3402" "3403" "3404" "3405" "3406" "3407" "3408" "3409" "3410"
## [3411] "3411" "3412" "3413" "3414" "3415" "3416" "3417" "3418" "3419" "3420"
## [3421] "3421" "3422" "3423" "3424" "3425" "3426" "3427" "3428" "3429" "3430"
## [3431] "3431" "3432" "3433" "3434" "3435" "3436" "3437" "3438" "3439" "3440"
## [3441] "3441" "3442" "3443" "3444" "3445" "3446" "3447" "3448" "3449" "3450"
## [3451] "3451" "3452" "3453" "3454" "3455" "3456" "3457" "3458" "3459" "3460"
## [3461] "3461" "3462" "3463" "3464" "3465" "3466" "3467" "3468" "3469" "3470"
## [3471] "3471" "3472" "3473" "3474" "3475" "3476" "3477" "3478" "3479" "3480"
## [3481] "3481" "3482" "3483" "3484" "3485" "3486" "3487" "3488" "3489" "3490"
## [3491] "3491" "3492" "3493" "3494" "3495" "3496" "3497" "3498" "3499" "3500"
## [3501] "3501" "3502" "3503" "3504" "3505" "3506" "3507" "3508" "3509" "3510"
## [3511] "3511" "3512" "3513" "3514" "3515" "3516" "3517" "3518" "3519" "3520"
## [3521] "3521" "3522" "3523" "3524" "3525" "3526" "3527" "3528" "3529" "3530"
## [3531] "3531" "3532" "3533" "3534" "3535" "3536" "3537" "3538" "3539" "3540"
## [3541] "3541" "3542" "3543" "3544" "3545" "3546" "3547" "3548" "3549" "3550"
## [3551] "3551" "3552" "3553" "3554" "3555" "3556" "3557" "3558" "3559" "3560"
## [3561] "3561" "3562" "3563" "3564" "3565" "3566" "3567" "3568" "3569" "3570"
## [3571] "3571" "3572" "3573" "3574" "3575" "3576" "3577" "3578" "3579" "3580"
## [3581] "3581" "3582" "3583" "3584" "3585" "3586" "3587" "3588" "3589" "3590"
## [3591] "3591" "3592" "3593" "3594" "3595" "3596" "3597" "3598" "3599" "3600"
## [3601] "3601" "3602" "3603" "3604" "3605" "3606" "3607" "3608" "3609" "3610"
## [3611] "3611" "3612" "3613" "3614" "3615" "3616" "3617" "3618" "3619" "3620"
## [3621] "3621" "3622" "3623" "3624" "3625" "3626" "3627" "3628" "3629" "3630"
## [3631] "3631" "3632" "3633" "3634" "3635" "3636" "3637" "3638" "3639" "3640"
## [3641] "3641" "3642" "3643" "3644" "3645" "3646" "3647" "3648" "3649" "3650"
## [3651] "3651" "3652" "3653" "3654" "3655" "3656" "3657" "3658" "3659" "3660"
## [3661] "3661" "3662" "3663" "3664" "3665" "3666" "3667" "3668" "3669" "3670"
## [3671] "3671" "3672" "3673" "3674" "3675" "3676" "3677" "3678" "3679" "3680"
## [3681] "3681" "3682" "3683" "3684" "3685" "3686" "3687" "3688" "3689" "3690"
## [3691] "3691" "3692" "3693" "3694" "3695" "3696" "3697" "3698" "3699" "3700"
## [3701] "3701" "3702" "3703" "3704" "3705" "3706" "3707" "3708" "3709" "3710"
## [3711] "3711" "3712" "3713" "3714" "3715" "3716" "3717" "3718" "3719" "3720"
## [3721] "3721" "3722" "3723" "3724" "3725" "3726" "3727" "3728" "3729" "3730"
## [3731] "3731" "3732" "3733" "3734" "3735" "3736" "3737" "3738" "3739" "3740"
## [3741] "3741" "3742" "3743" "3744" "3745" "3746" "3747" "3748" "3749" "3750"
## [3751] "3751" "3752" "3753" "3754" "3755" "3756" "3757" "3758" "3759" "3760"
## [3761] "3761" "3762" "3763" "3764" "3765" "3766" "3767" "3768" "3769" "3770"
## [3771] "3771" "3772" "3773" "3774" "3775" "3776" "3777" "3778" "3779" "3780"
## [3781] "3781" "3782" "3783" "3784" "3785" "3786" "3787" "3788" "3789" "3790"
## [3791] "3791" "3792" "3793" "3794" "3795" "3796" "3797" "3798" "3799" "3800"
## [3801] "3801" "3802" "3803" "3804" "3805" "3806" "3807" "3808" "3809" "3810"
## [3811] "3811" "3812" "3813" "3814" "3815" "3816" "3817" "3818" "3819" "3820"
## [3821] "3821" "3822" "3823" "3824" "3825" "3826" "3827" "3828" "3829" "3830"
## [3831] "3831" "3832" "3833" "3834" "3835" "3836" "3837" "3838" "3839" "3840"
## [3841] "3841" "3842" "3843" "3844" "3845" "3846" "3847" "3848" "3849" "3850"
## [3851] "3851" "3852" "3853" "3854" "3855" "3856" "3857" "3858" "3859" "3860"
## [3861] "3861" "3862" "3863" "3864" "3865" "3866" "3867" "3868" "3869" "3870"
## [3871] "3871" "3872" "3873" "3874" "3875" "3876" "3877" "3878" "3879" "3880"
## [3881] "3881" "3882" "3883" "3884" "3885" "3886" "3887" "3888" "3889" "3890"
## [3891] "3891" "3892" "3893" "3894" "3895" "3896" "3897" "3898" "3899" "3900"
## [3901] "3901" "3902" "3903" "3904" "3905" "3906" "3907" "3908" "3909" "3910"
## [3911] "3911" "3912" "3913" "3914" "3915" "3916" "3917" "3918" "3919" "3920"
## [3921] "3921" "3922" "3923" "3924" "3925" "3926" "3927" "3928" "3929" "3930"
## [3931] "3931" "3932" "3933" "3934" "3935" "3936" "3937" "3938" "3939" "3940"
## [3941] "3941" "3942" "3943" "3944" "3945" "3946" "3947" "3948" "3949" "3950"
## [3951] "3951" "3952" "3953" "3954" "3955" "3956" "3957" "3958" "3959" "3960"
## [3961] "3961" "3962" "3963" "3964" "3965" "3966" "3967" "3968" "3969" "3970"
## [3971] "3971" "3972" "3973" "3974" "3975" "3976" "3977" "3978" "3979" "3980"
## [3981] "3981" "3982" "3983" "3984" "3985" "3986" "3987" "3988" "3989" "3990"
## [3991] "3991" "3992" "3993" "3994" "3995" "3996" "3997" "3998" "3999" "4000"
## [4001] "4001" "4002" "4003" "4004" "4005" "4006" "4007" "4008" "4009" "4010"
## [4011] "4011" "4012" "4013" "4014" "4015" "4016" "4017" "4018" "4019" "4020"
## [4021] "4021" "4022" "4023" "4024" "4025" "4026" "4027" "4028" "4029" "4030"
## [4031] "4031" "4032" "4033" "4034" "4035" "4036" "4037" "4038" "4039" "4040"
## [4041] "4041" "4042" "4043" "4044" "4045" "4046" "4047" "4048" "4049" "4050"
## [4051] "4051" "4052" "4053" "4054" "4055" "4056" "4057" "4058" "4059" "4060"
## [4061] "4061" "4062" "4063" "4064" "4065" "4066" "4067" "4068" "4069" "4070"
## [4071] "4071" "4072" "4073" "4074" "4075" "4076" "4077" "4078" "4079" "4080"
## [4081] "4081" "4082" "4083" "4084" "4085" "4086" "4087" "4088" "4089" "4090"
## [4091] "4091" "4092" "4093" "4094" "4095" "4096" "4097" "4098" "4099" "4100"
## [4101] "4101" "4102" "4103" "4104" "4105" "4106" "4107" "4108" "4109" "4110"
## [4111] "4111" "4112" "4113" "4114" "4115" "4116" "4117" "4118" "4119" "4120"
## [4121] "4121" "4122" "4123" "4124" "4125" "4126" "4127" "4128" "4129" "4130"
## [4131] "4131" "4132" "4133" "4134" "4135" "4136" "4137" "4138" "4139" "4140"
## [4141] "4141" "4142" "4143" "4144" "4145" "4146" "4147" "4148" "4149" "4150"
## [4151] "4151" "4152" "4153" "4154" "4155" "4156" "4157" "4158" "4159" "4160"
## [4161] "4161" "4162" "4163" "4164" "4165" "4166" "4167" "4168" "4169" "4170"
## [4171] "4171" "4172" "4173" "4174" "4175" "4176" "4177" "4178" "4179" "4180"
## [4181] "4181" "4182" "4183" "4184" "4185" "4186" "4187" "4188" "4189" "4190"
## [4191] "4191" "4192" "4193" "4194" "4195" "4196" "4197" "4198" "4199" "4200"
## [4201] "4201" "4202" "4203" "4204" "4205" "4206" "4207" "4208" "4209" "4210"
## [4211] "4211" "4212" "4213" "4214" "4215" "4216" "4217" "4218" "4219" "4220"
## [4221] "4221" "4222" "4223" "4224" "4225" "4226" "4227" "4228" "4229" "4230"
## [4231] "4231" "4232" "4233" "4234" "4235" "4236" "4237" "4238" "4239" "4240"
## [4241] "4241" "4242" "4243" "4244" "4245" "4246" "4247" "4248" "4249" "4250"
## [4251] "4251" "4252" "4253" "4254" "4255" "4256" "4257" "4258" "4259" "4260"
## [4261] "4261" "4262" "4263" "4264" "4265" "4266" "4267" "4268" "4269" "4270"
## [4271] "4271" "4272" "4273" "4274" "4275" "4276" "4277" "4278" "4279" "4280"
## [4281] "4281" "4282" "4283" "4284" "4285" "4286" "4287" "4288" "4289" "4290"
## [4291] "4291" "4292" "4293" "4294" "4295" "4296" "4297" "4298" "4299" "4300"
## [4301] "4301" "4302" "4303" "4304" "4305" "4306" "4307" "4308" "4309" "4310"
## [4311] "4311" "4312" "4313" "4314" "4315" "4316" "4317" "4318" "4319" "4320"
## [4321] "4321" "4322" "4323" "4324" "4325" "4326" "4327" "4328" "4329" "4330"
## [4331] "4331" "4332" "4333" "4334" "4335" "4336" "4337" "4338" "4339" "4340"
## [4341] "4341" "4342" "4343" "4344" "4345" "4346" "4347" "4348" "4349" "4350"
## [4351] "4351" "4352" "4353" "4354" "4355" "4356" "4357" "4358" "4359" "4360"
## [4361] "4361" "4362" "4363" "4364" "4365" "4366" "4367" "4368" "4369" "4370"
## [4371] "4371" "4372" "4373" "4374" "4375" "4376" "4377" "4378" "4379" "4380"
## [4381] "4381" "4382" "4383" "4384" "4385" "4386" "4387" "4388" "4389" "4390"
## [4391] "4391" "4392" "4393" "4394" "4395" "4396" "4397" "4398" "4399" "4400"
## [4401] "4401" "4402" "4403" "4404" "4405" "4406" "4407" "4408" "4409" "4410"
## [4411] "4411" "4412" "4413" "4414" "4415" "4416" "4417" "4418" "4419" "4420"
## [4421] "4421" "4422" "4423" "4424" "4425" "4426" "4427" "4428" "4429" "4430"
## [4431] "4431" "4432" "4433" "4434" "4435" "4436" "4437" "4438" "4439" "4440"
## [4441] "4441" "4442" "4443" "4444" "4445" "4446" "4447" "4448" "4449" "4450"
## [4451] "4451" "4452" "4453" "4454" "4455" "4456" "4457" "4458" "4459" "4460"
## [4461] "4461" "4462" "4463" "4464" "4465" "4466" "4467" "4468" "4469" "4470"
## [4471] "4471" "4472" "4473" "4474" "4475" "4476" "4477" "4478" "4479" "4480"
## [4481] "4481" "4482" "4483" "4484" "4485" "4486" "4487" "4488" "4489" "4490"
## [4491] "4491" "4492" "4493" "4494" "4495" "4496" "4497" "4498" "4499" "4500"
## [4501] "4501" "4502" "4503" "4504" "4505" "4506" "4507" "4508" "4509" "4510"
## [4511] "4511" "4512" "4513" "4514" "4515" "4516" "4517" "4518" "4519" "4520"
## [4521] "4521" "4522" "4523" "4524" "4525" "4526" "4527" "4528" "4529" "4530"
## [4531] "4531" "4532" "4533" "4534" "4535" "4536" "4537" "4538" "4539" "4540"
## [4541] "4541" "4542" "4543" "4544" "4545" "4546" "4547" "4548" "4549" "4550"
## [4551] "4551" "4552" "4553" "4554" "4555" "4556" "4557" "4558" "4559" "4560"
## [4561] "4561" "4562" "4563" "4564" "4565" "4566" "4567" "4568" "4569" "4570"
## [4571] "4571" "4572" "4573" "4574" "4575" "4576" "4577" "4578" "4579" "4580"
## [4581] "4581" "4582" "4583" "4584" "4585" "4586" "4587" "4588" "4589" "4590"
## [4591] "4591" "4592" "4593" "4594" "4595" "4596" "4597" "4598" "4599" "4600"
## [4601] "4601" "4602" "4603" "4604" "4605" "4606" "4607" "4608" "4609" "4610"
## [4611] "4611" "4612" "4613" "4614" "4615" "4616" "4617" "4618" "4619" "4620"
## [4621] "4621" "4622" "4623" "4624" "4625" "4626" "4627" "4628" "4629" "4630"
## [4631] "4631" "4632" "4633" "4634" "4635" "4636" "4637" "4638" "4639" "4640"
## [4641] "4641" "4642" "4643" "4644" "4645" "4646" "4647" "4648" "4649" "4650"
## [4651] "4651" "4652" "4653" "4654" "4655" "4656" "4657" "4658" "4659" "4660"
## [4661] "4661" "4662" "4663" "4664" "4665" "4666" "4667" "4668" "4669" "4670"
## [4671] "4671" "4672" "4673" "4674" "4675" "4676" "4677" "4678" "4679" "4680"
## [4681] "4681" "4682" "4683" "4684" "4685" "4686" "4687" "4688" "4689" "4690"
## [4691] "4691" "4692" "4693" "4694" "4695" "4696" "4697" "4698" "4699" "4700"
## [4701] "4701" "4702" "4703" "4704" "4705" "4706" "4707" "4708" "4709" "4710"
## [4711] "4711" "4712" "4713" "4714" "4715" "4716" "4717" "4718" "4719" "4720"
## [4721] "4721" "4722" "4723" "4724" "4725" "4726" "4727" "4728" "4729" "4730"
## [4731] "4731" "4732" "4733" "4734" "4735" "4736" "4737" "4738" "4739" "4740"
## [4741] "4741" "4742" "4743" "4744" "4745" "4746" "4747" "4748" "4749" "4750"
## [4751] "4751" "4752" "4753" "4754" "4755" "4756" "4757" "4758" "4759" "4760"
## [4761] "4761" "4762" "4763" "4764" "4765" "4766" "4767" "4768" "4769" "4770"
## [4771] "4771" "4772" "4773" "4774" "4775" "4776" "4777" "4778" "4779" "4780"
## [4781] "4781" "4782" "4783" "4784" "4785" "4786" "4787" "4788" "4789" "4790"
## [4791] "4791" "4792" "4793" "4794" "4795" "4796" "4797" "4798" "4799" "4800"
## [4801] "4801" "4802" "4803" "4804" "4805" "4806" "4807" "4808" "4809" "4810"
## [4811] "4811" "4812" "4813" "4814" "4815" "4816" "4817" "4818" "4819" "4820"
## [4821] "4821" "4822" "4823" "4824" "4825" "4826" "4827" "4828" "4829" "4830"
## [4831] "4831" "4832" "4833" "4834" "4835" "4836" "4837" "4838" "4839" "4840"
## [4841] "4841" "4842" "4843" "4844" "4845" "4846" "4847" "4848" "4849" "4850"
## [4851] "4851" "4852" "4853" "4854" "4855" "4856" "4857" "4858" "4859" "4860"
## [4861] "4861" "4862" "4863" "4864" "4865" "4866" "4867" "4868" "4869" "4870"
## [4871] "4871" "4872" "4873" "4874" "4875" "4876" "4877" "4878" "4879" "4880"
## [4881] "4881" "4882" "4883" "4884" "4885" "4886" "4887" "4888" "4889" "4890"
## [4891] "4891" "4892" "4893" "4894" "4895" "4896" "4897" "4898" "4899" "4900"
## [4901] "4901" "4902" "4903" "4904" "4905" "4906" "4907" "4908" "4909" "4910"
## [4911] "4911" "4912" "4913" "4914" "4915" "4916" "4917" "4918" "4919" "4920"
## [4921] "4921" "4922" "4923" "4924" "4925" "4926" "4927" "4928" "4929" "4930"
## [4931] "4931" "4932" "4933" "4934" "4935" "4936" "4937" "4938" "4939" "4940"
## [4941] "4941" "4942" "4943" "4944" "4945" "4946" "4947" "4948" "4949" "4950"
## [4951] "4951" "4952" "4953" "4954" "4955" "4956" "4957" "4958" "4959" "4960"
## [4961] "4961" "4962" "4963" "4964" "4965" "4966" "4967" "4968" "4969" "4970"
## [4971] "4971" "4972" "4973" "4974" "4975" "4976" "4977" "4978" "4979" "4980"
## [4981] "4981" "4982" "4983" "4984" "4985" "4986" "4987" "4988" "4989" "4990"
## [4991] "4991" "4992" "4993" "4994" "4995" "4996" "4997" "4998" "4999" "5000"
## [5001] "5001" "5002" "5003" "5004" "5005" "5006" "5007" "5008" "5009" "5010"
## [5011] "5011" "5012" "5013" "5014" "5015" "5016" "5017" "5018" "5019" "5020"
## [5021] "5021" "5022" "5023" "5024" "5025" "5026" "5027" "5028" "5029" "5030"
## [5031] "5031" "5032" "5033" "5034" "5035" "5036" "5037" "5038" "5039" "5040"
## [5041] "5041" "5042" "5043" "5044" "5045" "5046" "5047" "5048" "5049" "5050"
## [5051] "5051" "5052" "5053" "5054" "5055" "5056" "5057" "5058" "5059" "5060"
## [5061] "5061" "5062" "5063" "5064" "5065" "5066" "5067" "5068" "5069" "5070"
## [5071] "5071" "5072" "5073" "5074" "5075" "5076" "5077" "5078" "5079" "5080"
## [5081] "5081" "5082" "5083" "5084" "5085" "5086" "5087" "5088" "5089" "5090"
## [5091] "5091" "5092" "5093" "5094" "5095" "5096" "5097" "5098" "5099" "5100"
## [5101] "5101" "5102" "5103" "5104" "5105" "5106" "5107" "5108" "5109" "5110"
## [5111] "5111" "5112" "5113" "5114" "5115" "5116" "5117" "5118" "5119" "5120"
## [5121] "5121" "5122" "5123" "5124" "5125" "5126" "5127" "5128" "5129" "5130"
## [5131] "5131" "5132" "5133" "5134" "5135" "5136" "5137" "5138" "5139" "5140"
## [5141] "5141" "5142" "5143" "5144" "5145" "5146" "5147" "5148" "5149" "5150"
## [5151] "5151" "5152" "5153" "5154" "5155" "5156" "5157" "5158" "5159" "5160"
## [5161] "5161" "5162" "5163" "5164" "5165" "5166" "5167" "5168" "5169" "5170"
## [5171] "5171" "5172" "5173" "5174" "5175" "5176" "5177" "5178" "5179" "5180"
## [5181] "5181" "5182" "5183" "5184" "5185" "5186" "5187" "5188" "5189" "5190"
## [5191] "5191" "5192" "5193" "5194" "5195" "5196" "5197" "5198" "5199" "5200"
## [5201] "5201" "5202" "5203" "5204" "5205" "5206" "5207" "5208" "5209" "5210"
## [5211] "5211" "5212" "5213" "5214" "5215" "5216" "5217" "5218" "5219" "5220"
## [5221] "5221" "5222" "5223" "5224" "5225" "5226" "5227" "5228" "5229" "5230"
## [5231] "5231" "5232" "5233" "5234" "5235" "5236" "5237" "5238" "5239" "5240"
## [5241] "5241" "5242" "5243" "5244" "5245" "5246" "5247" "5248" "5249" "5250"
## [5251] "5251" "5252" "5253" "5254" "5255" "5256" "5257" "5258" "5259" "5260"
## [5261] "5261" "5262" "5263" "5264" "5265" "5266" "5267" "5268" "5269" "5270"
## [5271] "5271" "5272" "5273" "5274" "5275" "5276" "5277" "5278" "5279" "5280"
## [5281] "5281" "5282" "5283" "5284" "5285" "5286" "5287" "5288" "5289" "5290"
## [5291] "5291" "5292" "5293" "5294" "5295" "5296" "5297" "5298" "5299" "5300"
## [5301] "5301" "5302" "5303" "5304" "5305" "5306" "5307" "5308" "5309" "5310"
## [5311] "5311" "5312" "5313" "5314" "5315" "5316" "5317" "5318" "5319" "5320"
## [5321] "5321" "5322" "5323" "5324" "5325" "5326" "5327" "5328" "5329" "5330"
## [5331] "5331" "5332" "5333" "5334" "5335" "5336" "5337" "5338" "5339" "5340"
## [5341] "5341" "5342" "5343" "5344" "5345" "5346" "5347" "5348" "5349" "5350"
## [5351] "5351" "5352" "5353" "5354" "5355" "5356" "5357" "5358" "5359" "5360"
## [5361] "5361" "5362" "5363" "5364" "5365" "5366" "5367" "5368" "5369" "5370"
## [5371] "5371" "5372" "5373" "5374" "5375" "5376" "5377" "5378" "5379" "5380"
## [5381] "5381" "5382" "5383" "5384" "5385" "5386" "5387" "5388" "5389" "5390"
## [5391] "5391" "5392" "5393" "5394" "5395" "5396" "5397" "5398" "5399" "5400"
## [5401] "5401" "5402" "5403" "5404" "5405" "5406" "5407" "5408" "5409" "5410"
## [5411] "5411" "5412" "5413" "5414" "5415" "5416" "5417" "5418" "5419" "5420"
## [5421] "5421" "5422" "5423" "5424" "5425" "5426" "5427" "5428" "5429" "5430"
## [5431] "5431" "5432" "5433" "5434" "5435" "5436" "5437" "5438" "5439" "5440"
## [5441] "5441" "5442" "5443" "5444" "5445" "5446" "5447" "5448" "5449" "5450"
## [5451] "5451" "5452" "5453" "5454" "5455" "5456" "5457" "5458" "5459" "5460"
## [5461] "5461" "5462" "5463" "5464" "5465" "5466" "5467" "5468" "5469" "5470"
## [5471] "5471" "5472" "5473" "5474" "5475" "5476" "5477" "5478" "5479" "5480"
## [5481] "5481" "5482" "5483" "5484" "5485" "5486" "5487" "5488" "5489" "5490"
## [5491] "5491" "5492" "5493" "5494" "5495" "5496" "5497" "5498" "5499" "5500"
## [5501] "5501" "5502" "5503" "5504" "5505" "5506" "5507" "5508" "5509" "5510"
## [5511] "5511" "5512" "5513" "5514" "5515" "5516" "5517" "5518" "5519" "5520"
## [5521] "5521" "5522" "5523" "5524" "5525" "5526" "5527" "5528" "5529" "5530"
## [5531] "5531" "5532" "5533" "5534" "5535" "5536" "5537" "5538" "5539" "5540"
## [5541] "5541" "5542" "5543" "5544" "5545" "5546" "5547" "5548" "5549" "5550"
## [5551] "5551" "5552" "5553" "5554" "5555" "5556" "5557" "5558" "5559" "5560"
## [5561] "5561" "5562" "5563" "5564" "5565" "5566" "5567" "5568" "5569" "5570"
## [5571] "5571" "5572" "5573" "5574" "5575" "5576" "5577" "5578" "5579" "5580"
## [5581] "5581" "5582" "5583" "5584" "5585" "5586" "5587" "5588" "5589" "5590"
## [5591] "5591" "5592" "5593" "5594" "5595" "5596" "5597" "5598" "5599" "5600"
## [5601] "5601" "5602" "5603" "5604" "5605" "5606" "5607" "5608" "5609" "5610"
## [5611] "5611" "5612" "5613" "5614" "5615" "5616" "5617" "5618" "5619" "5620"
## [5621] "5621" "5622" "5623" "5624" "5625" "5626" "5627" "5628" "5629" "5630"
## [5631] "5631" "5632" "5633" "5634" "5635" "5636" "5637" "5638" "5639" "5640"
## [5641] "5641" "5642" "5643" "5644" "5645" "5646" "5647" "5648" "5649" "5650"
## [5651] "5651" "5652" "5653" "5654" "5655" "5656" "5657" "5658" "5659" "5660"
## [5661] "5661" "5662" "5663" "5664" "5665" "5666" "5667" "5668" "5669" "5670"
## [5671] "5671" "5672" "5673" "5674" "5675" "5676" "5677" "5678" "5679" "5680"
## [5681] "5681" "5682" "5683" "5684" "5685" "5686" "5687" "5688" "5689" "5690"
## [5691] "5691" "5692" "5693" "5694" "5695" "5696" "5697" "5698" "5699" "5700"
## [5701] "5701" "5702" "5703" "5704" "5705" "5706" "5707" "5708" "5709" "5710"
## [5711] "5711" "5712" "5713" "5714" "5715" "5716" "5717" "5718" "5719" "5720"
## [5721] "5721" "5722" "5723" "5724" "5725" "5726" "5727" "5728" "5729" "5730"
## [5731] "5731" "5732" "5733" "5734" "5735" "5736" "5737" "5738" "5739" "5740"
## [5741] "5741" "5742" "5743" "5744" "5745" "5746" "5747" "5748" "5749" "5750"
## [5751] "5751" "5752" "5753" "5754" "5755" "5756" "5757" "5758" "5759" "5760"
## [5761] "5761" "5762" "5763" "5764" "5765" "5766" "5767" "5768" "5769" "5770"
## [5771] "5771" "5772" "5773" "5774" "5775" "5776" "5777" "5778" "5779" "5780"
## [5781] "5781" "5782" "5783" "5784" "5785" "5786" "5787" "5788" "5789" "5790"
## [5791] "5791" "5792" "5793" "5794" "5795" "5796" "5797" "5798" "5799" "5800"
## [5801] "5801" "5802" "5803" "5804" "5805" "5806" "5807" "5808" "5809" "5810"
## [5811] "5811" "5812" "5813" "5814" "5815" "5816" "5817" "5818" "5819" "5820"
## [5821] "5821" "5822" "5823" "5824" "5825" "5826" "5827" "5828" "5829" "5830"
## [5831] "5831" "5832" "5833" "5834" "5835" "5836" "5837" "5838" "5839" "5840"
## [5841] "5841" "5842" "5843" "5844" "5845" "5846" "5847" "5848" "5849" "5850"
## [5851] "5851" "5852" "5853" "5854" "5855" "5856" "5857" "5858" "5859" "5860"
## [5861] "5861" "5862" "5863" "5864" "5865" "5866" "5867" "5868" "5869" "5870"
## [5871] "5871" "5872" "5873" "5874" "5875" "5876" "5877" "5878" "5879" "5880"
## [5881] "5881" "5882" "5883" "5884" "5885" "5886" "5887" "5888" "5889" "5890"
## [5891] "5891" "5892" "5893" "5894" "5895" "5896" "5897" "5898" "5899" "5900"
## [5901] "5901" "5902" "5903" "5904" "5905" "5906" "5907" "5908" "5909" "5910"
## [5911] "5911" "5912" "5913" "5914" "5915" "5916" "5917" "5918" "5919" "5920"
## [5921] "5921" "5922" "5923" "5924" "5925" "5926" "5927" "5928" "5929" "5930"
## [5931] "5931" "5932" "5933" "5934" "5935" "5936" "5937" "5938" "5939" "5940"
## [5941] "5941" "5942" "5943" "5944" "5945" "5946" "5947" "5948" "5949" "5950"
## [5951] "5951" "5952" "5953" "5954" "5955" "5956" "5957" "5958" "5959" "5960"
## [5961] "5961" "5962" "5963" "5964" "5965" "5966" "5967" "5968" "5969" "5970"
## [5971] "5971" "5972" "5973" "5974" "5975" "5976" "5977" "5978" "5979" "5980"
## [5981] "5981" "5982" "5983" "5984" "5985" "5986" "5987" "5988" "5989" "5990"
## [5991] "5991" "5992" "5993" "5994" "5995" "5996" "5997" "5998" "5999" "6000"
## 
## $chain
## [1] "1" "2" "3" "4"
## 
## $variable
##  [1] "b_Intercept"            "b_sigma_Intercept"      "b_sFat"                
##  [4] "b_sAge"                 "b_breedckcs"            "b_breedpug"            
##  [7] "b_breedret"             "b_breedother"           "b_sexMale"             
## [10] "b_comorbidityyes"       "b_sigma_comorbidityyes" "sd_visit__Intercept"   
## [13] "alpha"                  "Intercept"              "Intercept_sigma"       
## [16] "r_visit[v0,Intercept]"  "r_visit[v1,Intercept]"  "prior_Intercept"       
## [19] "prior_b_sFat"           "prior_b_sAge"           "prior_b_breedckcs"     
## [22] "prior_b_breedpug"       "prior_b_breedret"       "prior_b_breedother"    
## [25] "prior_b_sexMale"        "prior_b_comorbidityyes" "prior_Intercept_sigma" 
## [28] "prior_alpha"            "prior_sd_visit"         "lprior"                
## [31] "lp__"                   "z_1[1,1]"               "z_1[1,2]"

ii create violin plot

mcmc_violin(posterior, pars = "b_sFat")

iii mcmc_pairs plot using bayesplot

color_scheme_set("viridis")
mcmc_trace(posterior, pars = "b_sFat") +
  
     theme_classic() +
     labs(title = "Trace plot for beta of standardised fat", 
         y = "b_sFat", 
         x = "Rank") +
         theme(axis.title.y = element_text(size=12, face="bold"), 
               axis.title.x = element_text(size=12, face="bold"),
               title = element_text(size=12, face="bold"),
               plot.title = element_text(hjust = 0.5),
               axis.text.x = element_text(color = "grey50", size = 12),
               axis.text.y = element_text(color = "grey8",size = 12))

mcmc_plot(model, type = 'rank_overlay', variable = "b_sFat") +

   theme_classic() +
     labs(title = "Trace rank plot for beta of standardised fat", 
         y = "b_sFat", 
         x = "Rank") +
         ylim(250, 350) +
         theme(axis.title.y = element_text(size=12, face="bold"), 
               axis.title.x = element_text(size=12, face="bold"),
               title = element_text(size=12, face="bold"),
               plot.title = element_text(hjust = 0.5),
               axis.text.x = element_text(color = "grey50", size = 12),
               axis.text.y = element_text(color = "grey8",size = 12))

## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.

8. Calculate 97% HPDI for fat

Usually better than the compatoability intervals given in the summary

draws <- as.matrix(model)
HPDI(draws[,3], 0.97) # 1st column is draws for bf

##      |0.97      0.97| 
## 0.01884829 0.58929953

9. Calculate R2 for model

bayes_R2(model, probs = c(0.015, 0.5, 0.985)) # 0.015, 0.5, 0.985 are the quantiles

##     Estimate  Est.Error       Q1.5       Q50     Q98.5
## R2 0.1835293 0.04804434 0.07906287 0.1837382 0.2893523

loo_R2(model, probs = c(0.015, 0.5, 0.985)) # 0.015, 0.5, 0.985 are the quantiles

## Warning: Some Pareto k diagnostic values are too high. See help('pareto-k-diagnostic') for details.

##      Estimate  Est.Error       Q1.5        Q50      Q98.5
## R2 -0.1560929 0.09813997 -0.3939027 -0.1521871 0.04130273

CHECKS ON MODEL

1. Basic check of simulations based on posterior distribution, versus the real data distribution

checks whether actual data is similar to simulated data. ### a. standard

color_scheme_set("blue")
pp_check(model, ndraws = 100) +

   theme_classic() +
  labs(title = "Posterior predictions of body fat model (skew-normal)", 
         y = "Density", 
         x = "Log hair cortisol (standardised") +
         theme(axis.title.y = element_text(size=12, face="bold"), 
               axis.title.x = element_text(size=12, face="bold"),
               title = element_text(size=12, face="bold"),
               plot.title = element_text(hjust = 0.5),
               axis.text.x = element_text(color = "grey50", size = 12),
               axis.text.y = element_text(color = "grey8",size = 12), ,
               legend.position = "inside", legend.position.inside = c(0.91, 0.83))

a. altered axes for comparison

color_scheme_set("blue")
pp_check(model, ndraws = 100) +

   theme_classic() +
  labs(title = "Posterior predictions of body fat model (skew-normal)", 
         y = "Density", 
         x = "Log hair cortisol (standardised") +
  xlim(-5, 7) +  
         theme(axis.title.y = element_text(size=12, face="bold"), 
               axis.title.x = element_text(size=12, face="bold"),
               title = element_text(size=12, face="bold"),
               plot.title = element_text(hjust = 0.5),
               axis.text.x = element_text(color = "grey50", size = 12),
               axis.text.y = element_text(color = "grey8",size = 12), ,
               legend.position = "inside", legend.position.inside = c(0.91, 0.83))

2. Check some individual draws versus observed using pp_check

par(mfrow = c(1,1))
pp_check(model, type = "hist", ndraws = 11, binwidth = 0.25) # separate histograms of 11 MCMC draws vs actual data

3. Other pp_check graphs

a. all

pp_check(model, type = "error_hist", ndraws = 11) # separate histograms of errors for 11 draws

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

pp_check(model, type = "scatter_avg", ndraws = 100) # scatter plot

pp_check(model, type = "stat_2d") #  scatterplot of joint posteriors

## Using all posterior draws for ppc type 'stat_2d' by default.
## Note: in most cases the default test statistic 'mean' is too weak to detect anything of interest.

# PPC functions for predictive checks based on (approximate) leave-one-out (LOO) cross-validation
pp_check(model, type = "loo_pit_overlay", ndraws = 1000) 

## Warning: Found 5 observations with a pareto_k > 0.67 in model '.x1'. We
## recommend to run more iterations to get at least about 2200 posterior draws to
## improve LOO-CV approximation accuracy.

## NOTE: The kernel density estimate assumes continuous observations and is not optimal for discrete observations.

b. for figure

i.

# separate histograms of 11 MCMC draws vs actual data
pp_check(model, type = "hist", ndraws = 11, binwidth = 0.25)  +

   theme_classic() +
     labs(title = "Actual data vs 11 MCMC draws for body fat model") +

         theme(axis.title.y = element_text(size=12, face="bold"), 
               axis.title.x = element_text(size=12, face="bold"),
               title = element_text(size=12, face="bold"),
               plot.title = element_text(hjust = 0.5),
               axis.text.x = element_text(color = "grey50", size = 12),
               axis.text.y = element_text(color = "grey8",size = 12))

ii.

#  scatterplot of joint posteriors
pp_check(model, type = "stat_2d")  +

   theme_classic() +
     labs(title = "Scatterplot of joint posteriors for body fat model") +

         theme(axis.title.y = element_text(size=12, face="bold"), 
               axis.title.x = element_text(size=12, face="bold"),
               title = element_text(size=12, face="bold"),
               plot.title = element_text(hjust = 0.5),
               axis.text.x = element_text(color = "grey50", size = 12),
               axis.text.y = element_text(color = "grey8",size = 12), ,
               legend.position = "inside", legend.position.inside = c(0.91, 0.13))

## Using all posterior draws for ppc type 'stat_2d' by default.

## Note: in most cases the default test statistic 'mean' is too weak to detect anything of interest.

pp_check(model, type = "loo_pit_overlay", ndraws = 1000) +

   theme_classic() +
     labs(title = "'LOO PIT overlay' plot for body fat model") +

         theme(axis.title.y = element_text(size=12, face="bold"), 
               axis.title.x = element_text(size=12, face="bold"),
               title = element_text(size=12, face="bold"),
               plot.title = element_text(hjust = 0.5),
               axis.text.x = element_text(color = "grey50", size = 12),
               axis.text.y = element_text(color = "grey8",size = 12),
               legend.position = "inside", legend.position.inside = c(0.91, 0.13))

## Warning: Found 7 observations with a pareto_k > 0.67 in model '.x1'. We
## recommend to run more iterations to get at least about 2200 posterior draws to
## improve LOO-CV approximation accuracy.

## NOTE: The kernel density estimate assumes continuous observations and is not optimal for discrete observations.

4. Error scatter average plot (recommended by Buerkner)

pp_check(model, type = "error_scatter_avg")

## Using all posterior draws for ppc type 'error_scatter_avg' by default.

5. Pairs plot

color_scheme_set("blue")
pairs(model, regex = TRUE)

color_scheme_set("pink")
mcmc_pairs(posterior,
           pars = c("Intercept", "Intercept_sigma",
                    "alpha", 
                    "b_sFat", "b_sAge",
                    "b_sexMale",
                    "b_comorbidityyes", "b_sigma_comorbidityyes"),
           off_diag_args = list(size = 0.5))

color_scheme_set("pink")
mcmc_pairs(posterior,
           pars = c("Intercept", "Intercept_sigma",
                    "alpha", 
                    "b_breedckcs", "b_breedpug",
                    "b_breedret", "b_breedother"),
                      off_diag_args = list(size = 0.5))

PSIS LOO-CV to check model performance

loo_model <- loo(model, moment_match = TRUE)
loo_model

## 
## Computed from 24000 by 55 log-likelihood matrix.
## 
##          Estimate   SE
## elpd_loo    -80.4  5.5
## p_loo        11.3  2.4
## looic       160.8 11.0
## ------
## MCSE of elpd_loo is 0.1.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.5, 1.2]).
## 
## All Pareto k estimates are good (k < 0.7).
## See help('pareto-k-diagnostic') for details.

AUTOMATED PRIOR SENSITIVITY USING THE PRIOR SENSE PACKAGE

1. Sensitivity check

First, check the sensitivity of the prior and likelihood to power-scaling. Posterior and posteriors resulting from power-scaling.

color_scheme_set("blue")
powerscale_sensitivity(model, variable = c("b_Intercept",
                                           "b_sFat",
                                           "b_sAge",
                                           "b_breedckcs", "b_breedother", "b_breedpug", "b_breedret",
                                           "b_sexMale",
                                           "b_comorbidityyes"))

## Sensitivity based on cjs_dist
## Prior selection: all priors
## Likelihood selection: all data
## 
##          variable prior likelihood diagnosis
##       b_Intercept 0.042      0.053         -
##            b_sFat 0.032      0.213         -
##            b_sAge 0.018      0.133         -
##       b_breedckcs 0.027      0.070         -
##      b_breedother 0.025      0.089         -
##        b_breedpug 0.025      0.112         -
##        b_breedret 0.027      0.092         -
##         b_sexMale 0.015      0.094         -
##  b_comorbidityyes 0.026      0.165         -

2. Now use bayestestR package to check priors are informative

check_prior(model, effects = "all")

##              Parameter Prior_Quality
## 1          b_Intercept   informative
## 2               b_sFat   informative
## 3               b_sAge   informative
## 4          b_breedckcs   informative
## 5           b_breedpug   informative
## 6           b_breedret   informative
## 7         b_breedother   informative
## 8            b_sexMale   informative
## 9     b_comorbidityyes   informative
## 10 sd_visit__Intercept   informative

CHECK PRIOR PREDICTION LINES FROM FINAL MODEL

1. Obtain draws of priors from final model

prior <- prior_draws(model)
prior %>% glimpse()

## Rows: 24,000
## Columns: 12
## $ Intercept        <dbl> 1.0290702041, 0.3306087896, 0.4839608205, 0.008573827…
## $ b_sFat           <dbl> 0.03607712, -0.54894592, 0.77995266, 0.83958415, -0.0…
## $ b_sAge           <dbl> -0.13587838, 0.21510924, 0.02527679, -0.81112013, 0.4…
## $ b_breedckcs      <dbl> 0.47432534, -2.09018501, 0.81528623, 0.42376328, -0.6…
## $ b_breedpug       <dbl> 1.79481210, -1.08212312, -2.04839099, 0.86380913, 0.5…
## $ b_breedret       <dbl> 0.9778925, -1.1427470, 0.9106284, -0.3588809, -1.0509…
## $ b_breedother     <dbl> -0.10398011, 1.11981942, 0.20578397, -2.64647193, 1.2…
## $ b_sexMale        <dbl> -1.260058872, 1.730863388, 1.910090289, -2.017788670,…
## $ b_comorbidityyes <dbl> -0.332723946, -0.002110186, 0.912710899, -0.232086474…
## $ Intercept_sigma  <dbl> -2.2490940, 3.9302854, -2.3171893, -5.9297446, 10.597…
## $ alpha            <dbl> 4.43100344, 3.73963395, 5.94892021, 5.50938705, 4.706…
## $ sd_visit         <dbl> 0.7810777, 0.1843761, 1.0139941, 1.1335228, 1.0756120…

2. Plot prior prediction lines for sAge

set.seed(5)

prior %>% 
  slice_sample(n = 50) %>% 
  rownames_to_column("draw") %>% 
  expand_grid(a = c(-2, 2)) %>% 
  mutate(c = Intercept + b_sAge * a) %>% 
  
  ggplot(aes(x = a, y = c)) +
  geom_line(aes(group = draw),
            color = "firebrick", alpha = .4) +
  labs(x = "Age (std)",
       y = "log(cort) (std)") +
  coord_cartesian(ylim = c(-4, 4)) +
  theme_bw() +
  theme(panel.grid = element_blank()) 

3. Plot prior prediction lines for sFat

a. basic

set.seed(5)

prior %>% 
  slice_sample(n = 50) %>% 
  rownames_to_column("draw") %>% 
  expand_grid(a = c(-2, 2)) %>% 
  mutate(c = Intercept + b_sFat * a) %>% 
  
  ggplot(aes(x = a, y = c)) +
  geom_line(aes(group = draw),
            color = "firebrick", alpha = .4) +
  labs(x = "Fat (std)",
       y = "log(cort) (std)") +
  coord_cartesian(ylim = c(-2, 2)) +
  theme_bw() +
  theme(panel.grid = element_blank()) 

b. for figure

i. prior prediction slopes

set.seed(5)

prior %>% 
  slice_sample(n = 50) %>% 
  rownames_to_column("draw") %>% 
  expand_grid(a = c(-2, 2)) %>% 
  mutate(c = Intercept + b_sFat * a) %>% 
  
  ggplot(aes(x = a, y = c)) +
  geom_line(aes(group = draw),
            color = "firebrick", alpha = .4) +
  labs(x = "Fat (standardised)",
       y = "Log hair cortisol (standardisedd)") +
  coord_cartesian(ylim = c(-2, 2)) +
     theme_classic() +
  theme(panel.grid = element_blank()) +

     labs(title = "Draws from prior for effect of body fat on hair cortisol") +  
         theme(axis.title.y = element_text(size=12, face="bold"), 
               axis.title.x = element_text(size=12, face="bold"),
               title = element_text(size=12, face="bold"),
               plot.title = element_text(hjust = 0.5),
               axis.text.x = element_text(color = "grey50", size = 12),
               axis.text.y = element_text(color = "grey8",size = 12))

i. posterior prediction slopes

set.seed(5)

post_draws <- draws
post_draws <- as.data.frame(post_draws)
post_draws %>% 
  slice_sample(n = 50) %>% 
  rownames_to_column("draw") %>% 
  expand_grid(a = c(-2, 2)) %>% 
  mutate(c = Intercept + b_sFat * a) %>% 
  
  ggplot(aes(x = a, y = c)) +
  geom_line(aes(group = draw),
            color = "steelblue3", alpha = .4) +
  labs(x = "Fat (standardised)",
       y = "Log hair cortisol (standardisedd)") +
  coord_cartesian(ylim = c(-2, 2)) +
     theme_classic() +
  theme(panel.grid = element_blank()) +

     labs(title = "Draws from posterior for effect of body fat on hair cortisol") +  
         theme(axis.title.y = element_text(size=12, face="bold"), 
               axis.title.x = element_text(size=12, face="bold"),
               title = element_text(size=12, face="bold"),
               plot.title = element_text(hjust = 0.5),
               axis.text.x = element_text(color = "grey50", size = 12),
               axis.text.y = element_text(color = "grey8",size = 12))

3. Plot prior prediction lines for comorbidity (yes)

set.seed(5)

prior %>% 
  slice_sample(n = 50) %>% 
  rownames_to_column("draw") %>% 
  expand_grid(a = c(0, 1)) %>% 
  mutate(c = Intercept + b_comorbidityyes * a) %>% 
  
  ggplot(aes(x = a, y = c)) +
  geom_line(aes(group = draw),
            color = "firebrick", alpha = .4) +
  geom_point(color = "firebrick", size = 2) +
  labs(x = "comorbidity",
       y = "log(cort) (std)") +
  coord_cartesian(ylim = c(-2, 2)) +
  theme_bw() +
  theme(panel.grid = element_blank()) 

4. Plot prior prediction lines for sex (male)

set.seed(5)

prior %>% 
  slice_sample(n = 50) %>% 
  rownames_to_column("draw") %>% 
  expand_grid(a = c(0, 1)) %>% 
  mutate(c = Intercept + b_sexMale * a) %>% 
  
  ggplot(aes(x = a, y = c)) +
  geom_line(aes(group = draw),
            color = "firebrick", alpha = .4) +
  geom_point(color = "firebrick", size = 2) +
  labs(x = "Sex breed",
       y = "log(cort) (std)") +
  coord_cartesian(ylim = c(-2, 2)) +
  theme_bw() +
  theme(panel.grid = element_blank()) 

5. Plot prior prediction lines for breedckcs

set.seed(5)

prior %>% 
  slice_sample(n = 50) %>% 
  rownames_to_column("draw") %>% 
  expand_grid(a = c(0, 1)) %>% 
  mutate(c = Intercept + b_breedckcs * a) %>% 
  
  ggplot(aes(x = a, y = c)) +
  geom_line(aes(group = draw),
            color = "firebrick", alpha = .4) +
  geom_point(color = "firebrick", size = 2) +
  labs(x = "CKCS breed",
       y = "log(cort) (std)") +
  coord_cartesian(ylim = c(-2, 2)) +
  theme_bw() +
  theme(panel.grid = element_blank()) 

6. Plot prior prediction lines for breedother

set.seed(5)

prior %>% 
  slice_sample(n = 50) %>% 
  rownames_to_column("draw") %>% 
  expand_grid(a = c(0, 1)) %>% 
  mutate(c = Intercept + b_breedother * a) %>% 
  
  ggplot(aes(x = a, y = c)) +
  geom_line(aes(group = draw),
            color = "firebrick", alpha = .4) +
  geom_point(color = "firebrick", size = 2) +
  labs(x = "Other breed",
       y = "log(cort) (std)") +
  coord_cartesian(ylim = c(-4, 4)) +
  theme_bw() +
  theme(panel.grid = element_blank()) 

7. Plot prior prediction lines for breedpug

set.seed(5)

prior %>% 
  slice_sample(n = 50) %>% 
  rownames_to_column("draw") %>% 
  expand_grid(a = c(0, 1)) %>% 
  mutate(c = Intercept + b_breedpug * a) %>% 
  
  ggplot(aes(x = a, y = c)) +
  geom_line(aes(group = draw),
            color = "firebrick", alpha = .4) +
  geom_point(color = "firebrick", size = 2) +
  labs(x = "Pug breed",
       y = "log(cort) (std)") +
  coord_cartesian(ylim = c(-4, 4)) +
  theme_bw() +
  theme(panel.grid = element_blank()) 

8. Plot prior prediction lines for breedret

set.seed(5)

prior %>% 
  slice_sample(n = 50) %>% 
  rownames_to_column("draw") %>% 
  expand_grid(a = c(0, 1)) %>% 
  mutate(c = Intercept + b_breedpug * a) %>% 
  
  ggplot(aes(x = a, y = c)) +
  geom_line(aes(group = draw),
            color = "firebrick", alpha = .4) +
  geom_point(color = "firebrick", size = 2) +
  labs(x = "Pug breed",
       y = "log(cort) (std)") +
  coord_cartesian(ylim = c(-4, 4)) +
  theme_bw() +
  theme(panel.grid = element_blank()) 

CHECK PRIOR PREDICTIVE DISTRIBUTION

1. Prior Predictive Distribution

Can simulate data just on the priors. Fit model but only consider prior when fitting model. If this looks reasonable, it helps to confirm that your priors were reasonable

set.seed(666)
model_priors_only <- brm(slgCort ~ sFat + sAge + breed + sex + comorbidity + (1 | visit),
                   family = skew_normal(),
                   prior = priors,
                   data = df2,
                   sample_prior = "only",
                   save_pars = save_pars(all =TRUE))

## Compiling Stan program...

## Trying to compile a simple C file

## Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
## using C compiler: ‘Apple clang version 17.0.0 (clang-1700.0.13.5)’
## using SDK: ‘MacOSX15.5.sdk’
## clang -arch arm64 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG   -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/Rcpp/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/unsupported"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/src/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppParallel/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/rstan/include" -DEIGEN_NO_DEBUG  -DBOOST_DISABLE_ASSERTS  -DBOOST_PENDING_INTEGER_LOG2_HPP  -DSTAN_THREADS  -DUSE_STANC3 -DSTRICT_R_HEADERS  -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION  -D_HAS_AUTO_PTR_ETC=0  -include '/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp'  -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1   -I/opt/R/arm64/include    -fPIC  -falign-functions=64 -Wall -g -O2  -c foo.c -o foo.o
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22:
## In file included from /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/Eigen/Dense:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/Eigen/Core:19:
## /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:679:10: fatal error: 'cmath' file not found
##   679 | #include <cmath>
##       |          ^~~~~~~
## 1 error generated.
## make: *** [foo.o] Error 1

## Start sampling

## 
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 1).
## Chain 1: 
## Chain 1: Gradient evaluation took 7.6e-05 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.76 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1: 
## Chain 1: 
## Chain 1: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 1: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 1: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 1: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 1: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 1: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 1: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 1: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 1: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 1: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 1: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 1: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 1: 
## Chain 1:  Elapsed Time: 0.034 seconds (Warm-up)
## Chain 1:                0.033 seconds (Sampling)
## Chain 1:                0.067 seconds (Total)
## Chain 1: 
## 
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 2).
## Chain 2: 
## Chain 2: Gradient evaluation took 4e-06 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.04 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2: 
## Chain 2: 
## Chain 2: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 2: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 2: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 2: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 2: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 2: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 2: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 2: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 2: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 2: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 2: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 2: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 2: 
## Chain 2:  Elapsed Time: 0.035 seconds (Warm-up)
## Chain 2:                0.028 seconds (Sampling)
## Chain 2:                0.063 seconds (Total)
## Chain 2: 
## 
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 3).
## Chain 3: 
## Chain 3: Gradient evaluation took 5e-06 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.05 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3: 
## Chain 3: 
## Chain 3: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 3: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 3: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 3: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 3: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 3: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 3: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 3: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 3: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 3: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 3: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 3: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 3: 
## Chain 3:  Elapsed Time: 0.038 seconds (Warm-up)
## Chain 3:                0.033 seconds (Sampling)
## Chain 3:                0.071 seconds (Total)
## Chain 3: 
## 
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 4).
## Chain 4: 
## Chain 4: Gradient evaluation took 3e-06 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.03 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4: 
## Chain 4: 
## Chain 4: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 4: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 4: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 4: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 4: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 4: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 4: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 4: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 4: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 4: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 4: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 4: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 4: 
## Chain 4:  Elapsed Time: 0.037 seconds (Warm-up)
## Chain 4:                0.027 seconds (Sampling)
## Chain 4:                0.064 seconds (Total)
## Chain 4:

## Warning: There were 1 divergent transitions after warmup. See
## https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
## to find out why this is a problem and how to eliminate them.

## Warning: Examine the pairs() plot to diagnose sampling problems

2. Check predictions against priors

color_scheme_set("pink")
pp_check(model_priors_only, ndraws = 100) +
  xlim(-10, 20) +

   theme_classic() +
     labs(title = "Prior predictions of body fat model versus data", 
         y = "Density", 
         x = "Log hair cortisol (standardised") +
         theme(axis.title.y = element_text(size=12, face="bold"), 
               axis.title.x = element_text(size=12, face="bold"),
               title = element_text(size=12, face="bold"),
               plot.title = element_text(hjust = 0.5),
               axis.text.x = element_text(color = "grey50", size = 12),
               axis.text.y = element_text(color = "grey8",size = 12),
               legend.position = "inside", legend.position.inside = c(0.91, 0.83))

## Warning: Removed 74 rows containing non-finite outside the scale range
## (`stat_density()`).

MANUAL POSTERIOR PREDICTIVE DISTRIBUTION CHECKS

NB Uses posterior_predict ## 1. Posterior predictive distribution plots for a continuous predictor variable

par(mfrow = c(2,2))
# plot the observed data
plot(slgCort ~ sFat, main = "Observed", xlab = "Fat (std)", ylab = "log hair cortisol (std)",
     data = df2, ylim = c(-3, 3), xlim = c(-3,3), col = "steelblue3") # This is the observed data

# use posterior predict to simulate predictions
ppd <- posterior_predict(model) 


# Plot 3 simulations of the data
plot(ppd[50,] ~ df2$sFat, main = "Simulated", xlab = "Fat (std)", ylab = "log hair cortisol (std)",
     ylim = c(-3,3), xlim = c(-3,3), col = "firebrick") 
plot(ppd[51,] ~ df2$sFat, main = "Simulated", xlab = "Fat (std)", ylab = "log hair cortisol (std)",
     ylim = c(-3,3), xlim = c(-3,3), col = "firebrick") 
plot(ppd[52,] ~ df2$sFat, main = "Simulated", xlab = "Fat (std)", ylab = "log hair cortisol (std)",
     ylim = c(-3,3), xlim = c(-3,3), col = "firebrick")

2. Posterior predictive distribution plots for a categorical predictor variable

par(mfrow = c(2,2))
vioplot(slgCort ~ comorbidity, data = df2, main = "Observed", col = "steelblue")
vioplot(ppd[sample(seq(1, dim(ppd)[1]), 1),] ~ df2$comorbidity, main = "PPD", col = "firebrick")
vioplot(ppd[sample(seq(1, dim(ppd)[1]), 1),] ~ df2$comorbidity, main = "PPD", col = "firebrick")
vioplot(ppd[sample(seq(1, dim(ppd)[1]), 1),] ~ df2$comorbidity, main = "PPD", col = "firebrick")

3. Posterior predictive distribition plots for a categorical predictor variable

par(mfrow = c(2,2))
vioplot(slgCort ~ sex, data = df2, main = "Observed", col = "steelblue3")
vioplot(ppd[sample(seq(1, dim(ppd)[1]), 1),] ~ df2$sex, main = "PPD", col = "firebrick")
vioplot(ppd[sample(seq(1, dim(ppd)[1]), 1),] ~ df2$sex, main = "PPD", col = "firebrick")
vioplot(ppd[sample(seq(1, dim(ppd)[1]), 1),] ~ df2$sex, main = "PPD", col = "firebrick")

4. Posterior predictive distribition plots for a categorical predictor variable with small group size

par(mfrow = c(2,2))
stripchart(slgCort ~ breed, vertical = TRUE, method = "jitter",
           col = "steelblue3", data = df2, pch = 20, main = "Observed")
stripchart(ppd[sample(seq(1, dim(ppd)[1]), 1),] ~ breed, vertical = TRUE, method = "jitter", col = "firebrick", data = df2, pch = 20, main = "PPD")
stripchart(ppd[sample(seq(1, dim(ppd)[1]), 1),] ~ breed, vertical = TRUE, method = "jitter", col = "firebrick", data = df2, pch = 20, main = "PPD")
stripchart(ppd[sample(seq(1, dim(ppd)[1]), 1),] ~ breed, vertical = TRUE, method = "jitter", col = "firebrick", data = df2, pch = 20, main = "PPD")

ANALYSING THE POSTERIOR DISTRIBUTION

1. Plot conditional effects from model

plot(conditional_effects(model), ask = FALSE)

b. format and plot conditional effects for weight_status

ce <- conditional_effects(model, "sFat", method = "posterior_predict", prob = 0.95)
plot(ce, plot = FALSE)[[1]] +
         theme_bw() +
         labs(title = "Conditional effect of body fat on hair cortisol") +
         labs(y = paste0("Log hair cortisol (standardised)")) +
         labs(x = paste0("Body fat (standardised)")) +
         theme(axis.title.y = element_text(size=12, face="bold"), 
               axis.title.x = element_text(size=12, face="bold"),
               title = element_text(size=12, face="bold"),
               plot.title = element_text(hjust = 0.5),
               axis.text.x = element_text(color = "grey50", size = 12),
               axis.text.y = element_text(color = "grey50", size = 12),
               legend.position = "inside", legend.position.inside = c(0.93, 0.80))

b. format and plot conditional effects for weight_status and comorbidity

ce <- conditional_effects(model, "sFat:comorbidity")
plot(ce, plot = FALSE)[[1]] +

           theme_bw() +
         labs(title = "Conditional effect of body fat on hair cortisol") +
         labs(y = paste0("Log hair cortisol (standardised)")) +
         labs(x = paste0("Body fat (standardised)")) +
         theme(axis.title.y = element_text(size=12, face="bold"), 
               axis.title.x = element_text(size=12, face="bold"),
               title = element_text(size=12, face="bold"),
               plot.title = element_text(hjust = 0.5),
               axis.text.x = element_text(color = "grey50", size = 12),
               axis.text.y = element_text(color = "grey50", size = 12),
               legend.position = "inside", legend.position.inside = c(0.91, 0.13))

mcmc_plot of model

color_scheme_set("blue")
mcmc_plot(model)

body fat vs prior

1. distributional

mcmc_plot(model,
          variable = c("b_sFat", "prior_b_sFat"))

2. density

mcmc_plot(model,
          variable = c("b_sFat", "prior_b_sFat"),
          type = "areas") +

   theme_classic() +
    labs(title = "Prior vs posterior distribution for body fat effect") +
         labs(y = "") +
         labs(x = paste0("Possible parameter values")) +
    scale_y_discrete(labels=c("prior_b_sFat" = "Prior", "b_sFat" = "Posterior"),
                     limits = c("prior_b_sFat", "b_sFat")) +
         theme(axis.title.y = element_text(size=12, face="bold"), 
               axis.title.x = element_text(size=12, face="bold"),
               title = element_text(size=12, face="bold"),
               plot.title = element_text(hjust = 0.5),
               axis.text.x = element_text(color = "grey50", size = 12),
               axis.text.y = element_text(color = "grey8",size = 12))

## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.

a.just parameters of beta variables

mcmc_plot(model,
          variable = "b_sexMale")

b. all parameters except

mcmc_plot(model, 
          variable = c("b_Intercept", "alpha",
         "b_sFat",
         "b_sAge",
         "b_breedckcs",
         "b_breedpug",
         "b_breedret",
         "b_breedother",
         "b_sexMale",
         "b_comorbidityyes"))

3. Plot all posterior distributions

posterior <- as.matrix(model)
mcmc_areas(posterior,
pars = c("b_Intercept", "alpha",
         "b_sFat",
         "b_sAge",
         "b_breedckcs",
         "b_breedpug",
         "b_breedret",
         "b_breedother",
         "b_sexMale",
         "b_comorbidityyes"),
# arbitrary threshold for shading probability mass
prob = 0.75)

4. plot posterior distribution for betas only

posterior <- as.matrix(model)
mcmc_areas(posterior,
    pars = c("b_sFat",
         "b_sAge",
         "b_breedckcs",
         "b_breedpug",
         "b_breedret",
         "b_breedother",
         "b_sexMale",
         "b_comorbidityyes"),
    prob = 0.75) # arbitrary threshold for shading probability mass

Plot posterior distribution for body fat percentage

posterior <- as.matrix(model)
mcmc_areas(posterior,
pars = "b_sFat",
# arbitrary threshold for shading probability mass
prob = 0.97) +
  
   theme_classic() +
     labs(title = "Posterior distribution for body fat effect", 
         y = "Density distribution", 
         x = "Possible parameter values") +
    scale_y_discrete(labels=("b_sFat" = "Beta for fat mass")) +
         theme(axis.title.y = element_text(size=12, face="bold"), 
               axis.title.x = element_text(size=12, face="bold"),
               title = element_text(size=12, face="bold"),
               plot.title = element_text(hjust = 0.5),
               axis.text.x = element_text(color = "grey50", size = 12),
               axis.text.y = element_text(color = "grey8",size = 12))

## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.

5. Describe the posterior visually

# Focus on describing posterior
hdi_range <- hdi(model, ci = c(0.65, 0.70, 0.80, 0.89, 0.95))
plot(hdi_range, show_intercept = T)

just sFat

# Focus on describing posterior
hdi_range <- hdi(model, ci = c(0.65, 0.70, 0.80, 0.89, 0.97), parameters = "sFat")
plot(hdi_range, show_intercept = T) +

    labs(title = "Posterior distribution for body fat effect") +
         labs(y = "Density distribution") +
         labs(x = "Possible parameter values") +
           theme(axis.title.y = element_text(size=12, face="bold"), 
               axis.title.x = element_text(size=12, face="bold"),
               title = element_text(size=12, face="bold"),
               plot.title = element_text(hjust = 0.5),
               axis.text.x = element_text(color = "grey50", size = 12),
               axis.text.y = element_text(color = "grey8",size = 12))

6. Effect of body fat on Log cortisol

##. a. Using ‘fitted’ which is a wrapper for posterior_epred (only returns error from model)

# Create new data for an average dog without a comorbidity
nd <- tibble(sFat = seq(from = -3.2, to = 3.2, length.out = 30), visit = "v0", breed = "mix", sex = "Female", comorbidity = "no", sAge = 0)

# now use `fitted()` to get the model-implied trajectories
pr1 <- fitted(model,
       newdata = nd) %>% 
  data.frame() %>% 
  bind_cols(nd) 

# Create new data for an average dog with a comorbidity
nd2 <- tibble(sFat = seq(from = -3.2, to = 3.2, length.out = 30), visit = "v0", breed = "mix", sex = "Female", comorbidity = "yes", sAge = 0)

# now use `fitted()` to get the model-implied trajectories
pr2 <- fitted(model,
       newdata = nd2) %>% 
  data.frame() %>% 
  bind_cols(nd2)


  # plot predictions for no comorbidity
  ggplot(data = pr1, aes(x = sFat)) +
  geom_smooth(aes(y = Estimate, ymin = Q2.5, ymax = Q97.5),
              stat = "identity",
              fill = "#F8766D", color = "#F8766D", alpha = 1/5, linewidth = 1.25) +
  
  # plot predictions for comorbidity
  geom_smooth(data = pr2, aes(y = Estimate, ymin = Q2.5, ymax = Q97.5),
              stat = "identity",
              fill = "#00BFC4", color = "#00BFC4", alpha = 1/5, linewidth = 1.25) + 
 
  
   # plot actual data
     geom_point(data = df2, 
             aes(y = slgCort, colour = comorbidity), 
             size = 3) +
 
   theme_bw() +
  labs(title = "Predicted effect of body fat on log hair cortisol",
       x = "Body fat (standardised)",
       y = "Log hair cortisol (standardised)") +
  theme(axis.title.y = element_text(size=14, face="bold"), 
        axis.title.x = element_text(size=14, face="bold"),
        title = element_text(size=14, face="bold"),
        plot.title = element_text(hjust = 0.5),
        legend.position = "inside", legend.position.inside = c(0.87, 0.13)) +
  coord_cartesian(ylim = c(-3, 3)) +
  scale_x_continuous(breaks=seq(-3, 3 , 1))

##. b. Using ‘predict’ which is a wrapper for posterior_pred (returns error from model and residual error)

# Create new data for an average dog without a comorbidity
nd <- tibble(sFat = seq(from = -3.2, to = 3.2, length.out = 30), visit = "v0", breed = "mix", sex = "Female", comorbidity = "no", sAge = 0)

# now use `fitted()` to get the model-implied trajectories
pr1 <- predict(model,
       newdata = nd) %>% 
  data.frame() %>% 
  bind_cols(nd) 

# Create new data for an average dog with a comorbidity
nd2 <- tibble(sFat = seq(from = -3.2, to = 3.2, length.out = 30), visit = "v0", breed = "mix", sex = "Female", comorbidity = "yes", sAge = 0)

# now use `fitted()` to get the model-implied trajectories
pr2 <- predict(model,
       newdata = nd2) %>% 
  data.frame() %>% 
  bind_cols(nd2)


  # plot predictions for no comorbidity
  ggplot(data = pr1, aes(x = sFat)) +
  geom_smooth(aes(y = Estimate, ymin = Q2.5, ymax = Q97.5),
              stat = "identity",
              fill = "#F8766D", color = "#F8766D", alpha = 1/5, linewidth = 1.25) +
  
  # plot predictions for comorbidity
  geom_smooth(data = pr2, aes(y = Estimate, ymin = Q2.5, ymax = Q97.5),
              stat = "identity",
              fill = "#00BFC4", color = "#00BFC4", alpha = 1/5, linewidth = 1.25) + 
 
  
   # plot actual data
     geom_point(data = df2, 
             aes(y = slgCort, colour = comorbidity), 
             size = 3) +
 
   theme_bw() +
  labs(title = "Predicted effect of body fat on log hair cortisol",
       x = "Body fat (standardised)",
       y = "Log hair cortisol (standardised)") +
  theme(axis.title.y = element_text(size=14, face="bold"), 
        axis.title.x = element_text(size=14, face="bold"),
        title = element_text(size=14, face="bold"),
        plot.title = element_text(hjust = 0.5),
        legend.position = "inside", legend.position.inside = c(0.87, 0.13)) +
  coord_cartesian(ylim = c(-3, 4)) +
  scale_x_continuous(breaks=seq(-3, 3 , 1))

HYPOTHESIS TESTS

1. Hypothesis test to check if mean association between cortisol and body fat (from draws) is >0

draws <- as.matrix(model)
mean(draws[,3] >0)

## [1] 0.9865417

2. Check 97% credible interval of with HPDI for body fat from draws

HPDI(draws[,3], prob=0.97)

##      |0.97      0.97| 
## 0.01884829 0.58929953

3. Hypothesis test that the effect of body fat is positive or negative

hypothesis(model, "sFat >0")

## Hypothesis Tests for class b:
##   Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio Post.Prob Star
## 1 (sFat) > 0     0.33      0.13      0.1     0.52       73.3      0.99    *
## ---
## 'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
## '*': For one-sided hypotheses, the posterior probability exceeds 95%;
## for two-sided hypotheses, the value tested against lies outside the 95%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.

hypothesis(model, "sFat <0")

## Hypothesis Tests for class b:
##   Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio Post.Prob Star
## 1 (sFat) < 0     0.33      0.13      0.1     0.52       0.01      0.01     
## ---
## 'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
## '*': For one-sided hypotheses, the posterior probability exceeds 95%;
## for two-sided hypotheses, the value tested against lies outside the 95%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.

POSTERIOR PREDICTIONS FOR DATA FROM NEW DATA

1. Plot the predicted effect of body fat mass on log hair cortisol

set.seed(666)
nd1 <- tibble(sFat = seq(-3, 3, by = 1), sAge = 0, visit = 0,
              sex = "Female", comorbidity = "no",
              breed = "mix",
              case_number = c("dog1", "dog2", "dog3","dog4", "dog5", "dog6", "dog7"))

p1 <-
  predict(model,
          resp = "slgCort",
          allow_new_levels = TRUE,
          newdata = nd1) %>% 
  data.frame() %>% 
  bind_cols(nd1) %>% 
  
  ggplot(aes(x = sFat, y = Estimate, ymin = Q2.5, ymax = Q97.5)) +
  
  geom_linerange(aes(ymin = Estimate - Est.Error, ymax = Estimate + Est.Error),
                 linewidth = 1, color = "#F8766D", alpha = 3/5) +
  geom_point(size = 5, color = "#F8766D") +
  theme_bw() +
  labs(title = "Predicted effect of body fat on hair cortisol",
       x = "Body fat (standardised)",
       y = "Log hair cortisol (standardised)") +
  theme(axis.title.y = element_text(size=12, face="bold"), 
        axis.title.x = element_text(size=12, face="bold"),
        title = element_text(size=12, face="bold"),
        plot.title = element_text(hjust = 0.5)) +
  coord_cartesian(ylim = c(-2.5, 2.5)) +
  scale_x_continuous(breaks=seq(-3, 3 , 1))

plot(p1)

7. Plot the predicted effect of body fat on log hair cortisol when comorbidity varies

a. create data for groups of dogs with different comorbidity status and predict results

set.seed(666)

nd3 <- tibble(sFat = seq(-3, 3, by = 1), sAge = 0, visit = 0,
              sex = "Female", comorbidity = "no",
              breed = "mix",
              case_number = c("dog1", "dog2", "dog3","dog4", "dog5", "dog6", "dog7"))

nd4 <- tibble(sFat = seq(-3, 3, by = 1), sAge = 0, visit = 0,
              sex = "Female", comorbidity = "yes",
              breed = "mix",
              case_number = c("dog1", "dog2", "dog3","dog4", "dog5", "dog6", "dog7"))


# predict outcome for comorbidity = no
pr3 <-
  predict(model,
          resp = "slgCort",
          allow_new_levels = TRUE,
          newdata = nd3) %>% 
  data.frame() %>% 
  bind_cols(nd3)

# predict outcome for comorbidity = yes
pr4 <-
  predict(model,
          resp = "slgCort",
          allow_new_levels = TRUE,
          newdata = nd4) %>% 
  data.frame() %>% 
  bind_cols(nd4)

# bind data for graphing
pr34 <-rbind(pr3, pr4)
pr34 <-  as.data.frame(pr34)

b. Plot the predicted estmates with estimated error

p3 <-  ggplot(pr34, aes(x = sFat, y = Estimate, , colour = comorbidity)) +
  
  geom_pointrange(aes(ymin = Estimate - Est.Error, ymax = Estimate + Est.Error),
                 linewidth = 1, size = 0.5, position = position_dodge2(width = 0.5, preserve = "single")) +
    theme_bw() +
    labs(title = "Predicted effect of body fat & comorbidity on hair cortisol") +
         labs(y = paste0("Predicted log hair cortisol (std)")) +
         labs(x = paste0("Body fat (std)")) +
         theme(axis.title.y = element_text(size=12, face="bold"), 
               axis.title.x = element_text(size=12, face="bold"),
               title = element_text(size=12, face="bold"),
               plot.title = element_text(hjust = 0.5)) +
          coord_cartesian(ylim = c(-2.5, 2.5), xlim = c(-3.5, 3.5)) +
         scale_x_continuous(limits = c(-3.5, 3.5), n.breaks = 7)

plot(p3)

4. Visualising the posterior of a model using numerical and graphical methods

a. basic (one dog only)

# create new dataframe which contains results of the first dog
new_data <- rbind(df2[1,], df2[1,])
# Now change one category to be different
new_data$sFat <- c(-1, 1)
# Visualise df to make sure it has worked
new_data

##   number   group visit season breed_group coat_colour  sex age comorbidity
## 1     c1 stopped    v0 winter         ret        dark Male  43         yes
## 2     c1 stopped    v0 winter         ret        dark Male  43         yes
##   fat_percent cortisol   lgCort breed sFat      sAge   slgCort
## 1    52.21393  4.92422 1.594166   ret   -1 -1.485464 0.3415375
## 2    52.21393  4.92422 1.594166   ret    1 -1.485464 0.3415375

# Now get predictions from the draws of the model
pred_new <- posterior_predict(model, newdata = new_data)


# Compare difference in means between the two categories
difference <- pred_new[,2] - pred_new[,1]


# Examine mean of difference
mean(difference)

## [1] 0.6415713

# View histogram of this
dens(difference)

# Create HPDI
HPDI(difference, 0.97)

##     |0.97     0.97| 
## -3.096552  4.522566

ii. plot effects

# Create data frame for the plot
pred_new <- data.frame(pred_new)
pred_new$draws <- seq(1:nrow(pred_new))
colnames(pred_new) <- c("sFat-1", "sFat+1", "draws")
dens_l <- melt(pred_new, id.vars = "draws", variable_name = "standardised_fat")


# Plot distribution of difference 
ggplot(dens_l, aes(value, fill = standardised_fat))  +
  geom_density(alpha = 1/2)  +
  xlim(-5, 10) +
    theme_bw() +
    labs(title = "Predicted log hair cortisol by body fat mass") +
         labs(y = paste0("Probability density")) +
         labs(x = paste0("Predicted Log Hair Cortisol (standardised)")) +
         theme(axis.title.y = element_text(size=12, face="bold"), 
               axis.title.x = element_text(size=12, face="bold"),
               title = element_text(size=12, face="bold"),
               plot.title = element_text(hjust = 0.5),
              legend.position = "inside", legend.position.inside = c(0.87, 0.87))

# create new dataframe which contains results of the first dog
new_data2 <- rbind(df2[1,], df2[1,])
# Now change one category to be different
new_data2$sFat <- c(-2, 2)
# Visualise df to make sure it has worked
new_data2

##   number   group visit season breed_group coat_colour  sex age comorbidity
## 1     c1 stopped    v0 winter         ret        dark Male  43         yes
## 2     c1 stopped    v0 winter         ret        dark Male  43         yes
##   fat_percent cortisol   lgCort breed sFat      sAge   slgCort
## 1    52.21393  4.92422 1.594166   ret   -2 -1.485464 0.3415375
## 2    52.21393  4.92422 1.594166   ret    2 -1.485464 0.3415375

# Now get predictions from the draws of the model
pred_new2 <- posterior_predict(model, newdata = new_data2)

# Create data frame for the plot
pred_new2 <- data.frame(pred_new2)
pred_new2$draws <- seq(1:nrow(pred_new2))
colnames(pred_new2) <- c("sFat-2", "sFat+2", "draws")
dens_l2 <- melt(pred_new2, id.vars = "draws", variable_name = "standardised_fat")


# Plot distribution of difference 
ggplot(dens_l2, aes(value, fill = standardised_fat))  +
  geom_density(alpha = 1/2)  +
  xlim(-6, 10) +
    theme_bw() +
    labs(title = "Predicted log hair cortisol by body fat mass") +
         labs(y = paste0("Probability density")) +
         labs(x = paste0("Predicted Log Hair Cortisol (standardised)")) +
         theme(axis.title.y = element_text(size=12, face="bold"), 
               axis.title.x = element_text(size=12, face="bold"),
               title = element_text(size=12, face="bold"),
               plot.title = element_text(hjust = 0.5),
              legend.position = "inside", legend.position.inside = c(0.87, 0.87))

b. Advanced… using all dogs in the model

# create new dataframe which contains results of all dogs
new_data1 <- df2
# Now change one category to be different
new_data1$sFat <- -1

# create new dataframe which contains results of the first dog
new_data2 <- df2
# Now change one category to be different
new_data2$sFat <- 1

# Now get predictions from the draws of the models
pred_nd1 <- posterior_predict(model, newdata = new_data1)
pred_nd2 <- posterior_predict(model, newdata = new_data2)
pred_diff <- pred_nd2 - pred_nd1
pred_diff <- data.frame(pred_diff)

# Create mean of differences for each column (dog) of the dataframe
pred_diff_mean <- apply(pred_diff, 2, mean)
# View histogram of mean differences
dens(pred_diff_mean)

# mean difference
mean(pred_diff_mean)

## [1] 0.6508525

# Create HPDI
HPDI(pred_diff_mean, 0.97)

##     |0.97     0.97| 
## 0.6335775 0.6681280

DO MODEL WITH GAUSSIAN DISTRIBUTION

1. Set priors

# Set individual priors
prior_int <- set_prior("normal(0, 0.5)", class = "Intercept")
prior_b <- set_prior("normal(0, 1)", class = "b")
prior_b_sex <- set_prior("normal(0, 1)", class = "b", coef = "sexMale")
prior_b_co <- set_prior("normal(0.25, 1)", class = "b", coef = "comorbidityyes")
prior_b_f <- set_prior("normal(0, 0.5)", class = "b", coef = "sFat")
prior_b_sAge <- set_prior("normal(0, 0.5)", class = "b", coef = "sAge")
prior_sd <- set_prior("normal(0, 1)", class = "sd")

# Combine priors into list
priors2 <- c(prior_int, prior_b, prior_b_sex, prior_b_co, prior_b_f, prior_b_sAge, prior_sd)

2. Run model2

Increased adapt_delta >0.8 (0.9999 here), as had divergent transitions

set.seed(666)
model2 <- brm(bf(slgCort ~ sFat + sAge + breed + sex + comorbidity + (1 | visit),
                   sigma ~ comorbidity),
                   family = gaussian(),
                   prior = priors2,
                   data = df2,
                   control=list(adapt_delta=0.9999, stepsize = 0.001, max_treedepth =15),
                   iter = 8000, warmup = 2000,
                   cores = 4,
                   save_pars = save_pars(all =TRUE),
                   sample_prior = TRUE)

## Compiling Stan program...

## Trying to compile a simple C file

## Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
## using C compiler: ‘Apple clang version 17.0.0 (clang-1700.0.13.5)’
## using SDK: ‘MacOSX15.5.sdk’
## clang -arch arm64 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG   -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/Rcpp/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/unsupported"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/src/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppParallel/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/rstan/include" -DEIGEN_NO_DEBUG  -DBOOST_DISABLE_ASSERTS  -DBOOST_PENDING_INTEGER_LOG2_HPP  -DSTAN_THREADS  -DUSE_STANC3 -DSTRICT_R_HEADERS  -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION  -D_HAS_AUTO_PTR_ETC=0  -include '/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp'  -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1   -I/opt/R/arm64/include    -fPIC  -falign-functions=64 -Wall -g -O2  -c foo.c -o foo.o
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22:
## In file included from /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/Eigen/Dense:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/Eigen/Core:19:
## /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:679:10: fatal error: 'cmath' file not found
##   679 | #include <cmath>
##       |          ^~~~~~~
## 1 error generated.
## make: *** [foo.o] Error 1

## Start sampling

## Found more than one class "stanfit" in cache; using the first, from namespace 'rethinking'

## Also defined by 'rstan'

## Found more than one class "stanfit" in cache; using the first, from namespace 'rethinking'

## Also defined by 'rstan'

## Found more than one class "stanfit" in cache; using the first, from namespace 'rethinking'

## Also defined by 'rstan'

## Found more than one class "stanfit" in cache; using the first, from namespace 'rethinking'

## Also defined by 'rstan'

3. Get summary of model2

summary(model2)

##  Family: gaussian 
##   Links: mu = identity; sigma = log 
## Formula: slgCort ~ sFat + sAge + breed + sex + comorbidity + (1 | visit) 
##          sigma ~ comorbidity
##    Data: df2 (Number of observations: 55) 
##   Draws: 4 chains, each with iter = 8000; warmup = 2000; thin = 1;
##          total post-warmup draws = 24000
## 
## Multilevel Hyperparameters:
## ~visit (Number of levels: 2) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.39      0.37     0.01     1.39 1.00     8922     9456
## 
## Regression Coefficients:
##                      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## Intercept               -0.80      0.45    -1.69     0.06 1.00    12410
## sigma_Intercept         -0.90      0.33    -1.47    -0.18 1.00    12297
## sFat                     0.25      0.13    -0.03     0.50 1.00    14898
## sAge                     0.06      0.15    -0.24     0.36 1.00    21911
## breedckcs               -0.06      0.50    -1.04     0.94 1.00    19989
## breedpug                 0.11      0.48    -0.83     1.07 1.00    14806
## breedret                -0.18      0.37    -0.90     0.55 1.00    14480
## breedother               0.39      0.39    -0.38     1.15 1.00    12968
## sexMale                  0.22      0.28    -0.33     0.76 1.00    20566
## comorbidityyes           0.84      0.26     0.32     1.34 1.00    19495
## sigma_comorbidityyes     1.04      0.36     0.28     1.67 1.00    13377
##                      Tail_ESS
## Intercept               15550
## sigma_Intercept         13064
## sFat                    14651
## sAge                    17194
## breedckcs               16647
## breedpug                17360
## breedret                16487
## breedother              16384
## sexMale                 16867
## comorbidityyes          16270
## sigma_comorbidityyes    14482
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

4. PSIS LOO-CV to check model performance

loo_model2 <- loo(model2, moment_match = TRUE)
loo_model2

## 
## Computed from 24000 by 55 log-likelihood matrix.
## 
##          Estimate   SE
## elpd_loo    -82.7  5.5
## p_loo        11.2  2.2
## looic       165.4 10.9
## ------
## MCSE of elpd_loo is 0.0.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.5, 1.2]).
## 
## All Pareto k estimates are good (k < 0.7).
## See help('pareto-k-diagnostic') for details.

5. Basic check of simulations based on posterior distribution, versus the real data distribution

checks whether actual data is similar to simulated data.

color_scheme_set("blue")
pp_check(model2, ndraws = 100) +

   theme_classic() +
   labs(title = "Posterior predictions of body fat model (normal)", 
         y = "Density", 
         x = "Log hair cortisol (standardised") +
    xlim(-5, 7) +
         theme(axis.title.y = element_text(size=12, face="bold"), 
               axis.title.x = element_text(size=12, face="bold"),
               title = element_text(size=12, face="bold"),
               plot.title = element_text(hjust = 0.5),
               axis.text.x = element_text(color = "grey50", size = 12),
               axis.text.y = element_text(color = "grey8",size = 12), ,
               legend.position = "inside", legend.position.inside = c(0.91, 0.83))

C. DO MODEL WITH Student’s T DISTRIBUTION

1. Set priors3

# Set individual priors
prior_int <- set_prior("normal(0, 0.5)", class = "Intercept")
prior_b <- set_prior("normal(0, 1)", class = "b")
prior_b_sex <- set_prior("normal(0, 1)", class = "b", coef = "sexMale")
prior_b_co <- set_prior("normal(0.25, 1)", class = "b", coef = "comorbidityyes")
prior_b_f <- set_prior("normal(0, 0.5)", class = "b", coef = "sFat")
prior_b_sAge <- set_prior("normal(0, 0.5)", class = "b", coef = "sAge")
prior_sd <- set_prior("normal(0, 1)", class = "sd")

# Combine priors into list
priors3 <- c(prior_int, prior_b, prior_b_sex, prior_b_co, prior_b_f, prior_b_sAge, prior_sd)

2. Run model3

Increased adapt_delta >0.8 (0.9999 here), as had divergent transitions

set.seed(666)
model3 <- brm(bf(slgCort ~ sFat + sAge + breed + sex + comorbidity + (1 | visit),
                   sigma ~ comorbidity),
                   family = student(),
                   prior = priors3,
                   data = df2,
                   control=list(adapt_delta=0.9999, stepsize = 0.001, max_treedepth =15),
                   iter = 8000, warmup = 2000,
                   cores = 4,
                   save_pars = save_pars(all =TRUE),
                   sample_prior = TRUE)

## Compiling Stan program...

## Trying to compile a simple C file

## Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
## using C compiler: ‘Apple clang version 17.0.0 (clang-1700.0.13.5)’
## using SDK: ‘MacOSX15.5.sdk’
## clang -arch arm64 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG   -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/Rcpp/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/unsupported"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/src/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppParallel/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/rstan/include" -DEIGEN_NO_DEBUG  -DBOOST_DISABLE_ASSERTS  -DBOOST_PENDING_INTEGER_LOG2_HPP  -DSTAN_THREADS  -DUSE_STANC3 -DSTRICT_R_HEADERS  -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION  -D_HAS_AUTO_PTR_ETC=0  -include '/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp'  -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1   -I/opt/R/arm64/include    -fPIC  -falign-functions=64 -Wall -g -O2  -c foo.c -o foo.o
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22:
## In file included from /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/Eigen/Dense:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/Eigen/Core:19:
## /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:679:10: fatal error: 'cmath' file not found
##   679 | #include <cmath>
##       |          ^~~~~~~
## 1 error generated.
## make: *** [foo.o] Error 1

## Start sampling

## Found more than one class "stanfit" in cache; using the first, from namespace 'rethinking'

## Also defined by 'rstan'

## Found more than one class "stanfit" in cache; using the first, from namespace 'rethinking'

## Also defined by 'rstan'

## Found more than one class "stanfit" in cache; using the first, from namespace 'rethinking'

## Also defined by 'rstan'

## Found more than one class "stanfit" in cache; using the first, from namespace 'rethinking'

## Also defined by 'rstan'

3. Get summary of model3

summary(model3)

##  Family: student 
##   Links: mu = identity; sigma = log; nu = identity 
## Formula: slgCort ~ sFat + sAge + breed + sex + comorbidity + (1 | visit) 
##          sigma ~ comorbidity
##    Data: df2 (Number of observations: 55) 
##   Draws: 4 chains, each with iter = 8000; warmup = 2000; thin = 1;
##          total post-warmup draws = 24000
## 
## Multilevel Hyperparameters:
## ~visit (Number of levels: 2) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.39      0.38     0.01     1.42 1.00    10092     9944
## 
## Regression Coefficients:
##                      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## Intercept               -0.79      0.45    -1.67     0.09 1.00    13581
## sigma_Intercept         -0.93      0.34    -1.53    -0.19 1.00    13282
## sFat                     0.24      0.14    -0.05     0.49 1.00    14641
## sAge                     0.04      0.16    -0.27     0.34 1.00    20236
## breedckcs               -0.09      0.50    -1.07     0.89 1.00    20540
## breedpug                 0.06      0.50    -0.90     1.08 1.00    15827
## breedret                -0.18      0.37    -0.90     0.55 1.00    16563
## breedother               0.38      0.40    -0.40     1.16 1.00    13904
## sexMale                  0.21      0.28    -0.34     0.75 1.00    21938
## comorbidityyes           0.82      0.27     0.28     1.33 1.00    17736
## sigma_comorbidityyes     1.01      0.37     0.23     1.69 1.00    14175
##                      Tail_ESS
## Intercept               15545
## sigma_Intercept         14262
## sFat                    14787
## sAge                    16248
## breedckcs               17568
## breedpug                16307
## breedret                17050
## breedother              16013
## sexMale                 17538
## comorbidityyes          15625
## sigma_comorbidityyes    15339
## 
## Further Distributional Parameters:
##    Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## nu    22.58     13.97     4.99    58.45 1.00    21794    14260
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

4.PSIS LOO-CV to check model performance

loo_model3 <- loo(model3, moment_match = TRUE)
loo_model3

## 
## Computed from 24000 by 55 log-likelihood matrix.
## 
##          Estimate   SE
## elpd_loo    -83.0  5.5
## p_loo        11.3  2.0
## looic       166.1 10.9
## ------
## MCSE of elpd_loo is 0.1.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.5, 1.1]).
## 
## All Pareto k estimates are good (k < 0.7).
## See help('pareto-k-diagnostic') for details.

5. Compare performance of models with different likelihhod functions

loo_compare(loo_model, loo_model2, loo_model3)

##        elpd_diff se_diff
## model   0.0       0.0   
## model2 -2.3       2.0   
## model3 -2.7       2.0

1. Basic check of simulations based on posterior distribution, versus the real data distribution

checks whether actual data is similar to simulated data.

color_scheme_set("blue")
pp_check(model3, ndraws = 100) +

   theme_classic() +
     labs(title = "Posterior predictions of body fat model (Student's t)", 
         y = "Density", 
         x = "Log hair cortisol (standardised") +
     xlim(-5, 7) +
         theme(axis.title.y = element_text(size=12, face="bold"), 
               axis.title.x = element_text(size=12, face="bold"),
               title = element_text(size=12, face="bold"),
               plot.title = element_text(hjust = 0.5),
               axis.text.x = element_text(color = "grey50", size = 12),
               axis.text.y = element_text(color = "grey8",size = 12), ,
               legend.position = "inside", legend.position.inside = c(0.91, 0.83))

C. DO MODEL WITH Student’s T DISTRIBUTION

1. Set priors

# Set individual priors
prior_int <- set_prior("normal(0, 0.5)", class = "Intercept")
prior_sig <- set_prior("exponential(1)", class = "sigma")
prior_b <- set_prior("normal(0, 1)", class = "b")
prior_b_sex <- set_prior("normal(0, 1)", class = "b", coef = "sexMale")
prior_b_co <- set_prior("normal(0.25, 1)", class = "b", coef = "comorbidityyes")
prior_b_f <- set_prior("normal(0, 0.5)", class = "b", coef = "sFat")
prior_b_sAge <- set_prior("normal(0, 0.5)", class = "b", coef = "sAge")
prior_sd <- set_prior("normal(0, 1)", class = "sd")
prior_alpha <- set_prior("normal(4, 2)", class = "alpha")

# Combine priors into list
priors4 <- c(prior_int, prior_sig, prior_b, prior_b_sex, prior_b_co, prior_b_f, prior_b_sAge, prior_sd, prior_alpha)

2. Run model4

Increased adapt_delta >0.8 (0.9999 here), as had divergent transitions

set.seed(666)
model4 <- brm(slgCort ~ sFat + sAge + breed + sex + comorbidity + (1 | visit),
                   family = skew_normal(),
                   prior = priors4,
                   data = df2,
                   control=list(adapt_delta=0.9999, stepsize = 0.001, max_treedepth =15),
                   iter = 8000, warmup = 2000,
                   cores = 4,
                   save_pars = save_pars(all =TRUE),
                   sample_prior = TRUE)

## Compiling Stan program...

## Trying to compile a simple C file

## Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
## using C compiler: ‘Apple clang version 17.0.0 (clang-1700.0.13.5)’
## using SDK: ‘MacOSX15.5.sdk’
## clang -arch arm64 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG   -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/Rcpp/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/unsupported"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/src/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppParallel/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/rstan/include" -DEIGEN_NO_DEBUG  -DBOOST_DISABLE_ASSERTS  -DBOOST_PENDING_INTEGER_LOG2_HPP  -DSTAN_THREADS  -DUSE_STANC3 -DSTRICT_R_HEADERS  -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION  -D_HAS_AUTO_PTR_ETC=0  -include '/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp'  -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1   -I/opt/R/arm64/include    -fPIC  -falign-functions=64 -Wall -g -O2  -c foo.c -o foo.o
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22:
## In file included from /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/Eigen/Dense:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/Eigen/Core:19:
## /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:679:10: fatal error: 'cmath' file not found
##   679 | #include <cmath>
##       |          ^~~~~~~
## 1 error generated.
## make: *** [foo.o] Error 1

## Start sampling

## Found more than one class "stanfit" in cache; using the first, from namespace 'rethinking'

## Also defined by 'rstan'

## Found more than one class "stanfit" in cache; using the first, from namespace 'rethinking'

## Also defined by 'rstan'

## Found more than one class "stanfit" in cache; using the first, from namespace 'rethinking'

## Also defined by 'rstan'

## Found more than one class "stanfit" in cache; using the first, from namespace 'rethinking'

## Also defined by 'rstan'

3. Get summary of model3

summary(model4)

##  Family: skew_normal 
##   Links: mu = identity; sigma = identity; alpha = identity 
## Formula: slgCort ~ sFat + sAge + breed + sex + comorbidity + (1 | visit) 
##    Data: df2 (Number of observations: 55) 
##   Draws: 4 chains, each with iter = 8000; warmup = 2000; thin = 1;
##          total post-warmup draws = 24000
## 
## Multilevel Hyperparameters:
## ~visit (Number of levels: 2) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.40      0.38     0.01     1.40 1.00     8219     9090
## 
## Regression Coefficients:
##                Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept         -0.27      0.50    -1.28     0.70 1.00    13870    14908
## sFat               0.18      0.17    -0.16     0.51 1.00    18969    17548
## sAge               0.01      0.15    -0.29     0.30 1.00    20284    16401
## breedckcs         -0.12      0.41    -0.94     0.68 1.00    17551    16049
## breedpug           0.19      0.56    -0.92     1.26 1.00    17030    15958
## breedret          -0.18      0.34    -0.83     0.49 1.00    15928    16747
## breedother         0.02      0.35    -0.65     0.73 1.00    14746    15376
## sexMale            0.15      0.25    -0.35     0.65 1.00    21009    16334
## comorbidityyes     0.40      0.36    -0.26     1.14 1.00    18046    15850
## 
## Further Distributional Parameters:
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     1.08      0.12     0.88     1.33 1.00    17460    16467
## alpha     3.93      1.63     1.04     7.47 1.00    14707    14861
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

4.PSIS LOO-CV to check model performance

loo_model4 <- loo(model4, moment_match = TRUE)
loo_model4

## 
## Computed from 24000 by 55 log-likelihood matrix.
## 
##          Estimate   SE
## elpd_loo    -85.1  5.7
## p_loo         9.7  2.3
## looic       170.1 11.3
## ------
## MCSE of elpd_loo is 0.1.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.6, 1.2]).
## 
## All Pareto k estimates are good (k < 0.7).
## See help('pareto-k-diagnostic') for details.

5. Compare performance of models with different likelihhod functions

loo_compare(loo_model, loo_model4)

##        elpd_diff se_diff
## model   0.0       0.0   
## model4 -4.7       2.2

6. Basic check of simulations based on posterior distribution, versus the real data distribution

checks whether actual data is similar to simulated data.

color_scheme_set("blue")
pp_check(model4, ndraws = 100) +

   theme_classic() +
     labs(title = "Posterior predictions of body fat model (Student's t)", 
         y = "Density", 
         x = "Log hair cortisol (standardised") +
         theme(axis.title.y = element_text(size=12, face="bold"), 
               axis.title.x = element_text(size=12, face="bold"),
               title = element_text(size=12, face="bold"),
               plot.title = element_text(hjust = 0.5),
               axis.text.x = element_text(color = "grey50", size = 12),
               axis.text.y = element_text(color = "grey8",size = 12), ,
               legend.position = "inside", legend.position.inside = c(0.91, 0.83))